🚀 go-pugleaf

BUG REPORT 17

BEHAVIOR:

fatal signal 4 produced upon compilation

COMPILER VERSION: GNU C++ compiler driver, version 1.25.0

INPUT FILES:

given below in shar file

COMMAND LINE: g++ -v -g rmdif.cc

FILE NAMES:

tm.h = tm-vax.h
md = vax.md

MACHINE: vax 11/780 with BRL UNIX (4.2 BSD)

SHAR FILE:
----------------------------------------------------------------------------
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	bug/Stream.h
#	bug/debug.h
#	bug/doubleCbuffer.h
#	bug/doubleComplex.h
#	bug/doubleFft.h
#	bug/lmath.h
#	bug/overload.h
#	bug/rmdif.cc
#	bug/std.h
# This archive created: Fri Jan  6 12:10:33 1989
export PATH; PATH=/bin:$PATH
if test -f 'bug/Stream.h'
then
	echo shar: will not over-write existing file "'bug/Stream.h'"
else
cat << \SHAR_EOF > 'bug/Stream.h'
//
// Stream.h : local augmented version of stream.h
//

#ifndef _Stream_h_
#define _Stream_h_

#include <overload.h>
#include <stream.h>

struct istream_init {
  istream_init()
    { set_File_error_handler(fatal_File_error_handler); }
} istream_init_dummy;

inline char* FileExt(const char* fname, const char* fext) {
  if ( !index(fname,'.') ) {
    char* c = new char[strlen(fname) + strlen(fext) + 2];
    strcpy(c,fname);
    strcat(c,".");
    strcat(c,fext);
    return c;
  }
  else
    return fname;
}

#endif
SHAR_EOF
fi # end of overwriting check
if test -f 'bug/debug.h'
then
	echo shar: will not over-write existing file "'bug/debug.h'"
else
cat << \SHAR_EOF > 'bug/debug.h'
//
// debug.h : code debugging aids
//

#ifndef _debug_h_
#define _debug_h_

#define show(x) cout << #x ## " = " << x << "\n"

#endif

SHAR_EOF
fi # end of overwriting check
if test -f 'bug/doubleCbuffer.h'
then
	echo shar: will not over-write existing file "'bug/doubleCbuffer.h'"
else
cat << \SHAR_EOF > 'bug/doubleCbuffer.h'
//
// doubleCbuffer.h : circular buffer class definition 
//

#ifndef _doubleCbuffer_h_
#define _doubleCbuffer_h_

#include <stream.h>
#include <error.h>

//
// doubleCbuffer : circular buffer class
//
//	supports
//	- run time buffer size specification
//	- buffer value load operator =
//	- first value extraction / assignment operator *
//	- Nth value extraction / assignment operator []
//	- buffer shift ++, --, +=, -=
//	- standard output operator <<
//
//	notes
//	- size 0 buffers are fully supported, except that
//	  requests for the contents of these buffers (via
//	  `*' and `[]') return garbage.
//	- subscripts of the `[]' operator are not checked
//	  for validity and will return garbage if illegal
//	  ( not 0 <= x < size ).
//

class doubleCbuffer {
  int	bsize;				// buffer size
  double *base,				// buffer base
    *end,				// end of buffer
    *start;				// current beginning

public:
  doubleCbuffer(int i) {
    bsize = i; 
    base = start = new double[bsize ? bsize : 1]; 
    end = base + bsize; }
  doubleCbuffer(doubleCbuffer& c)
    { doubleCbuffer(c.bsize,c); }
  doubleCbuffer(int i, doubleCbuffer& c) {
    bsize = i;
    base = start = new double[bsize ? bsize : 1]; 
    end = base + bsize;
    double *p = base; 
    for ( int j=0; j < c.bsize; j++ ) *p++ = c[j]; }
  doubleCbuffer(int i, double *q) {
    bsize = i;
    base = start = new double[bsize ? bsize : 1]; 
    end = base + bsize;
    for ( double *p=base; p < end; ) *p++ = *q++; }
  doubleCbuffer(int i, double d) {
    bsize = i;
    base = start = new double[bsize ? bsize : 1]; 
    end = base + bsize;
    for ( double *p=base; p < end; ) *p++ = d; }
  doubleCbuffer(int i, istream& s) {
    bsize = i;
    base = start = new double[bsize ? bsize : 1]; 
    end = base + bsize;
    for ( double *p=base; p < end; ) s >> *p++; 
    if ( !s ) warning("doubleCbuffer(istream&) : error on input"); }

  ~doubleCbuffer()
    { delete base; }

  doubleCbuffer& operator=(doubleCbuffer& c) {
    if ( this != &c ) {
      start = base;
      double *p=base;
      for ( int i=0; i < c.bsize; i++ ) *p++ = c[i];
    }
    return *this; }
  doubleCbuffer& operator=(double *q) {
    start = base;
    for ( double *p=base; p < end; ) *p++ = *q++;
    return *this; }

  double& operator*()
    { return *start; }
  double& operator[](int i)
    { return start+i >= end ? *(start+i-bsize) : *(start+i); }

  doubleCbuffer& operator++()
    { if ( ++start >= end ) start = base; return *this; }
  doubleCbuffer& operator--()
    { if ( --start < base ) start = end-1; return *this; }

  doubleCbuffer& operator+=(int i) {
    start += i; 
    if ( bsize )
      while ( start >= end ) start -= bsize; 
    return *this; }
  doubleCbuffer& operator-=(int i) {
    start -= i; 
    if ( bsize )
      while ( start < base ) start += bsize; 
    return *this; }

  friend int size(doubleCbuffer& c)
    { return c.bsize; }

  friend ostream& operator<<(ostream& s, doubleCbuffer& c) {
    s << "[";
    for ( int i=0; i < c.bsize; i++ ) s << " " << c[i]; 
    s << " ]";
    return s; }
};

#endif
SHAR_EOF
fi # end of overwriting check
if test -f 'bug/doubleComplex.h'
then
	echo shar: will not over-write existing file "'bug/doubleComplex.h'"
else
cat << \SHAR_EOF > 'bug/doubleComplex.h'
//
// doubleComplex.h : prototypable complex class
//

#ifndef _doubleComplex_h
#define _doubleComplex_h 1

#include <stream.h>
#include <std.h>

class doubleComplex
{
  double              re;
  double              im;

public:
                   doubleComplex();
                   doubleComplex(doubleComplex& c);
                   doubleComplex(double r, double i = 0.0);

                   ~doubleComplex();

  doubleComplex&      operator =  (doubleComplex& y);

  friend int       operator == (doubleComplex& x, doubleComplex& y);
  friend int       operator == (doubleComplex& x, double y);

  friend int       operator != (doubleComplex& x, doubleComplex& y);
  friend int       operator != (doubleComplex& x, double y);

  doubleComplex       operator +  ();
  doubleComplex       operator -  ();
  doubleComplex       operator ~  ();

  friend doubleComplex   operator +  (doubleComplex& x, doubleComplex& y);
  friend doubleComplex   operator -  (doubleComplex& x, doubleComplex& y);
  friend doubleComplex   operator *  (doubleComplex& x, doubleComplex& y);
  friend doubleComplex   operator /  (doubleComplex& x, doubleComplex& y);

  doubleComplex&         operator += (doubleComplex& y); 
  doubleComplex&         operator -= (doubleComplex& y); 
  doubleComplex&         operator *= (doubleComplex& y); 
  doubleComplex&         operator /= (doubleComplex& y); 

  friend double&   real(doubleComplex& x);
  friend double&   imag(doubleComplex& x);
  friend double    norm(doubleComplex& x);

  friend istream&  operator >> (istream& s, doubleComplex& x);
  friend ostream&  operator << (ostream& s, doubleComplex& x);

  void             error(char* msg);
};

// error handlers

extern  void default_doubleComplex_error_handler(char*);
extern  one_arg_error_handler_t doubleComplex_error_handler;

extern  one_arg_error_handler_t 
        set_doubleComplex_error_handler(one_arg_error_handler_t f);

//#ifdef __OPTIMIZE__

inline doubleComplex:: doubleComplex() {}
inline doubleComplex::~doubleComplex() {}

inline doubleComplex::doubleComplex(double r, double i = 0.0)
{
  re = r;
  im = i;
}

inline doubleComplex::doubleComplex(doubleComplex& x)
{
  re = x.re;
  im = x.im;
}

inline doubleComplex& doubleComplex::operator = (doubleComplex& x)
{
  re = x.re;
  im = x.im;
  return *this;
}

inline int operator == (doubleComplex& x, doubleComplex& y)
{
  return x.re == y.re && x.im == y.im;
}

inline int operator == (doubleComplex& x, double y)
{
  return x.im == 0.0 && x.re == y;
}

inline int operator != (doubleComplex& x, doubleComplex& y)
{
  return x.re != y.re || x.im != y.im;
}

inline int operator != (doubleComplex& x, double y)
{
  return x.im != 0.0 || x.re != y;
}

inline doubleComplex doubleComplex::operator + ()
{
  return *this;
}

inline doubleComplex doubleComplex::operator - ()
{
  return doubleComplex(-re, -im);
}

inline doubleComplex doubleComplex::operator ~ ()
{
  return doubleComplex(re, -im);
}

inline doubleComplex operator + (doubleComplex& x, doubleComplex& y)
{
  return doubleComplex(x.re + y.re, x.im + y.im);
}

inline doubleComplex operator - (doubleComplex& x, doubleComplex& y)
{
  return doubleComplex(x.re - y.re, x.im - y.im);
}

inline doubleComplex operator * (doubleComplex& x, doubleComplex& y)
{
  return doubleComplex(x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re);
}

inline doubleComplex& doubleComplex::operator += (doubleComplex& y)
{
  re += y.re;
  im += y.im;
  return *this;
}

inline doubleComplex& doubleComplex::operator -= (doubleComplex& y)
{
  re -= y.re;
  im -= y.im;
  return *this;
}

inline doubleComplex& doubleComplex::operator *= (doubleComplex& y)
{
  double r = re * y.re - im * y.im;
  im = re * y.im + im * y.re;
  re = r;
  return *this;
}

inline double& real(doubleComplex& x)
{
  return x.re;
}

inline double& imag(doubleComplex& x)
{
  return x.im;
}

inline double norm(doubleComplex& x)
{
  return (x.re * x.re + x.im * x.im);
}

//#endif

#endif
SHAR_EOF
fi # end of overwriting check
if test -f 'bug/doubleFft.h'
then
	echo shar: will not over-write existing file "'bug/doubleFft.h'"
else
cat << \SHAR_EOF > 'bug/doubleFft.h'
//
// doubleFft.h : double fft functions
//

#ifndef _doubleFft_h
#define _doubleFft_h

#include <doubleComplex.h>

overload fft;

// twiddle factor array
class doubleTwiddle {
  int fsize;				// fft size
  int dsize;				// data size: lengthof(data) 
  doubleComplex* data;			// twiddle factors

public:
  doubleTwiddle(int N, double = 1.0);		// constructor: size & stage gain
  ~doubleTwiddle() { delete data; }

  friend doubleComplex* fft(doubleComplex*, const doubleTwiddle&);
};

#endif
SHAR_EOF
fi # end of overwriting check
if test -f 'bug/lmath.h'
then
	echo shar: will not over-write existing file "'bug/lmath.h'"
else
cat << \SHAR_EOF > 'bug/lmath.h'

//
// lmath.h : local math definition file
//

#ifndef _lmath_h_
#define _lmath_h_

#include <overload.h>

inline int max(int a,int b)
  { return a > b ? a : b; }
inline int min(int a,int b)
  { return a < b ? a : b; }

inline long max(long a,long b)
  { return a > b ? a : b; }
inline long min(long a,long b)
  { return a < b ? a : b; }

inline double max(double a,double b)
  { return a > b ? a : b; }
inline double min(double a,double b)
  { return a < b ? a : b; }

#endif
SHAR_EOF
fi # end of overwriting check
if test -f 'bug/overload.h'
then
	echo shar: will not over-write existing file "'bug/overload.h'"
else
cat << \SHAR_EOF > 'bug/overload.h'
//
// overload.h : single repository for overload statements
//
//	this single repository avoids warning messages of multiple
//	overloads.
//

#ifndef _overload_h_
#define _overload_h_

overload warning;
overload error;
overload abort;

overload square;
overload min;
overload max;

#endif
SHAR_EOF
fi # end of overwriting check
if test -f 'bug/rmdif.cc'
then
	echo shar: will not over-write existing file "'bug/rmdif.cc'"
else
cat << \SHAR_EOF > 'bug/rmdif.cc'
//
// rmdif: rapid measurement of digital instantaneous frequency
//	  taken from Griffiths: ASSP-23, No. 2, pp. 207-22
//
// input:  sampled data file
// output: sample number / dectected signal pairs

// header inclusion
#include <Stream.h>
#include <math.h>
#include <lmath.h>
#include <debug.h>
#include <Random.h>
#include <ACG.h>

// type definition -- set type in following:
#include <doubleComplex.h>
#include <doubleFft.h>
#include <doubleCbuffer.h>

// type definition -- change the following in the rest of the code:
// double
// doubleComplex
// doubleTwiddle 
// doubleCbuffer

// external data
char usage[] = "\
usage: rmdif alpha filter_length snr signal_gain output_interval \n\
             fft_length detection_threshold";

// external definitions
#define LMAX	255

// macro definitions
#define argument(n, scanstring, x)				\
	if ( argc > n && sscanf(argv[n], scanstring, &x) != 1 )	\
	  error(usage)

// function declarations
double noise(double x, double signal_gain, double noise_gain);
void lms(double x, double mu, double* g, doubleCbuffer z);
void mspec(int M, int L, int N, double* g, doubleTwiddle W, double threshold);

// main function
main(int argc, char** argv) {
  int L = 4,				// filter length
    I = 4,				// output interval
    N = 256;				// fft length

  double alpha = 1.0,			// lms parameter
    snr = 1000.0,			// signal-to-noise ratio
    sgain = 1.0,			// signal gain
    threshold = 50.;			// detection threshold

  if ( argc > 7 )
    error(usage);
  argument(1, "%f", alpha);
  argument(2, "%d", L);
  argument(3, "%f", snr);
  argument(4, "%f", sgain);
  argument(5, "%d", I);
  argument(6, "%d", N);

  if ( alpha <= 0. || L < 2 || L >= LMAX || sgain <= 0. || I <= 0 || N < 4 )
    error("rmdif: parameter out of bounds");

  double ngain 				// noise gain
      = (snr >= 1000.0 ? 0.0 : sgain * pow(10, -snr/20.) / sqrt(2.)),
    x;					// input data

  double mu				// lms paramter
      = alpha / max( L * ( sgain * sgain / 2. + ngain * ngain ), .0001 ),
    X,					// noisy input data
    T = 1. / threshold;			// spectrum detection threshold

  double g[L+1];			// filter coefficients
  g[0] = 1.;

  doubleCbuffer z(L,0.);		// filter delay line

  doubleTwiddle W(N,0.5);		// twiddle factors (stage gain = .5)

  int i = I,				// output count
    m = 1;				// sample count

  // process data
  for ( ; cin >> x; m++ ) {
    X = noise(x, sgain, ngain);
    lms(X, mu, g, z);
    if ( --i <= 0 ) {
      i = I;
      mspec(m, L+1, N, g, W, T);
    }
  }
}

double noise(double x, double signal_gain, double noise_gain) {
  static ACG gen(10,20);
  static Normal rnd(0,1);

  return signal_gain * x + noise_gain * rnd();
}

void lms(double x, double mu, double* g, doubleCbuffer z) {
  for ( int i=0; i < size(z); i++ )
    g[i+1] = g[i+1] + mu * ( x + g[i+1] * z[i] );
  *--z = x;
}

overload MGmin;
inline int MGmin(double T, double g1, double g2, double g3) {
  return g1 < T && g1 < g2 && g1 < g3;
}

inline int MGmin(double T, double g1, double g2) {
  return g1 < T && g1 < g2;
}

void mspec(int M, int L, int N, double* g, doubleTwiddle W, double T) {
  doubleComplex G[N];			// fft workspace

  // form fft of filter coefficients
  for ( int i=0; i < L; i++ )		// initialize fft workspace
    G[i] = g[i];
  for ( ; i < N; i++ )
    G[i] = 0;
  fft(G,W);				// perform fft operation

  double MG[N];				// set magnitude of fft
  for ( i=0; i < N; i++ )
    MG[i] = norm(G[i]);

  if ( MGmin(T, MG[0], MG[1]) )		// find all local minima
    cout << M << " " << 0.0 << "\n";
  for ( i=1; i < N; i++ )
    if ( MGmin(T, MG[i], MG[i-1], MG[i+1]) )
      cout << M << " " << double(i)/N << "\n";
  if ( MGmin(T, MG[N-1], MG[N-2]) )
    cout << M << " " << double(N/2-1)/N << "\n";
}
SHAR_EOF
fi # end of overwriting check
if test -f 'bug/std.h'
then
	echo shar: will not over-write existing file "'bug/std.h'"
else
cat << \SHAR_EOF > 'bug/std.h'
//
// std.h : local version of std.h
//

#include <overload.h>
#include "/usr/local/lib/g++-include/std.h"

SHAR_EOF
fi # end of overwriting check
#	End of shell archive
exit 0