pse.c

这是数字信号处理方面的一些源码

/*********************************************************************

PSE.C - 用FFT实现功率谱估计

程序输入一个加入噪声的正弦信号, 采用周期图法来实现功率谱分析

参数如下:
 slice = 每次FFT的长度,缺省为256;
 numav = FFT取平均的组数,缺省为16;
 ovlap = 每次FFT运算交叉的点数,缺省为128;
 wtype = 选用窗函数的类型: 1--5.
 ntype = 选用噪声信号类型, 0:无噪声; 1:均匀噪声; 2:高斯噪声

使用子程序为:
fft            基2 DIF FFT 子程序
ifft           基2 DIF IFFT 子程序
log2           基2求对数函数
addnoise       在输入信号上加入均匀分布或高斯分布的噪声
gaussian       产生均值为零,单位方差的高斯分布随机数
uniform        产生零均值,均匀分布的随机数(-0.5 到 0.5 之间)
ham            Hamming 窗函数
han            Hanning 窗函数
triang         triangle 窗函数
black          Blackman 窗函数
harris         4 term Blackman-Harris 窗函数
draw_image     绘图子程序
***********************************************************************/

#include    <math.h>
#include    <stdlib.h>
#include    <stdio.h>
#include    <string.h>
#include    <conio.h>
#include    <graphics.h>

#define MAX(a,b)    (((a) > (b)) ? (a) : (b))
#define    PI	    (4.0*atan(1.0))

/* COMPLEX STRUCTURE */
typedef struct {
    float real, imag;
} COMPLEX;

 void fft(COMPLEX *,int);
 void ifft(COMPLEX *,int);
 int  log2(unsigned int x);
 double *addnoise(double *Indata,unsigned length,int type,double nmult);
 double gaussian(void);
 double uniform(void);
 void ham(COMPLEX *x, int n);
 void han(COMPLEX *x, int n);
 void triang(COMPLEX *,int n);
 void black(COMPLEX *x, int n);
 void harris(COMPLEX *x, int n);
 void draw_image(double *x,int m,char *title1,char *title2,
		char *xdis1,char *xdis2,int dis_type);
/*********************************************************************/

int numav=16;     /* 每次FFT的长度*/
int slice=256;    /* FFT取平均的组数*/
int ovlap=128;    /* 每次FFT运算交叉的点数 */
int wtype=5;      /* 选用窗函数的类型: 0--5 */
int ntype=2;      /* 选用噪声信号类型, 0:无噪声; 1:均匀噪声; 2:高斯噪声 */


void main(void)
{
  int        i, j, m, insize,index;
  int        k,numests,estsize;
  double     a,*mag, *sigfloat;
  double     tempflt;
  char	     title[80],tmp[20];
  COMPLEX    *samp;

  estsize = (slice-ovlap)*(numav-1) + slice;
  printf("The size of each estimate is: %d\n", estsize);
  sigfloat = (double *) calloc(estsize,sizeof(double));
  mag = (double *) calloc(slice,sizeof(double));
  samp = (COMPLEX *) calloc(slice,sizeof(COMPLEX));
  if (!mag || !samp){
    printf("\n  Memory allocation error.\n");
    exit(1);
  }
  /* input signal data series */
  for (i=0; i<estsize; i++)
      sigfloat[i] = cos(6.0*1.0*PI*(double)i*2.0/slice);
  for (i=0; i<estsize; i++)
      sigfloat[i] = sigfloat[i]+sin(35.0*1.0*PI*(double)i*2.0/slice);

  if(ntype==1)  addnoise(sigfloat,estsize,0,0.2); /* 加入均匀噪声 */
  else if(ntype==2) addnoise(sigfloat,estsize,1,0.2); /* 加入高斯噪声 */

  strcpy(title,"The Input Signal Data (x[n])");
  draw_image(sigfloat,estsize,title,"The Magnitude value","0",
	     itoa(estsize,tmp,10),0);

  /* calculate the periodical signal estimation value */
  for (k=0; k<slice; k++) mag[k] = 0;
  for (j=0; j<numav; j++){
     for (k=0; k<slice; k++){
	index = j*(slice-ovlap) + k;
	samp[k].real = sigfloat[index];
	samp[k].imag = 0;
     }
     switch (wtype) {
	case 1: ham(samp,slice);
		break;
	case 2: han(samp,slice);
		break;
	case 3: triang(samp,slice);
		break;
	case 4: black(samp,slice);
		break;
	case 5: harris(samp,slice);
		break;
	default:
		break;
     }

     m = log2(slice);
     fft(samp,m);
     a = (double) slice*slice;
     for (k=0; k<slice; k++){
	tempflt  = samp[k].real * samp[k].real;
	tempflt += samp[k].imag * samp[k].imag;
	tempflt = tempflt/a;
	mag[k] += tempflt;
     }
     printf("*");
  }

/*  Take log after averaging the magnitudes.  */

  for (k=0; k<slice/2; k++){
     mag[k] = mag[k]/numav;
     mag[k] = 10*log10(MAX(mag[k],1.e-14));
  }

  strcpy(title,"The Estimation Reault (X[m])");
  draw_image(mag,slice/2,title,"The Magnitude value (dB)","0",
	     "PI",0);

  free(mag);
  free(sigfloat);
  free(samp);
}

/************************************************************************
fft - 基2 DIF FFT 子程序

输入参数:
	 COMPLEX *x : FFT 输入和输出数据区指针;
	      int m : FFT 长度 ( length = 2^m );
输出参数:
	 输出数据放在 x 所指的输入数据区.
	 无输出参数.

void fft(COMPLEX *x, int m)
*************************************************************************/

void fft(COMPLEX *x,int m)
{
    static COMPLEX *w;           /* used to store the w complex array */
    static int mstore = 0;       /* stores m for future reference */
    static int n = 1;            /* length of fft stored for future */

    COMPLEX u,temp,tm;
    COMPLEX *xi,*xip,*xj,*wptr;

    int i,j,k,l,le,windex;

    double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;

    if(m != mstore) {

/* free previously allocated storage and set new m */

	if(mstore != 0) free(w);
	mstore = m;
	if(m == 0) return;       /* if m=0 then done */

/* n = 2**m = fft length */

	n = 1 << m;
	le = n/2;

/* allocate the storage for w */

	w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
	if(!w) {
	    printf("\nUnable to allocate complex W array\n");
	    exit(1);
	}

/* calculate the w values recursively */

	arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
	wrecur_real = w_real = cos(arg);
	wrecur_imag = w_imag = -sin(arg);
	xj = w;
	for (j = 1 ; j < le ; j++) {
	    xj->real = (float)wrecur_real;
	    xj->imag = (float)wrecur_imag;
	    xj++;
	    wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
	    wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
	    wrecur_real = wtemp_real;
	}
    }

/* start fft */

    le = n;
    windex = 1;
    for (l = 0 ; l < m ; l++) {
	le = le/2;

/* first iteration with no multiplies */

	for(i = 0 ; i < n ; i = i + 2*le) {
	    xi = x + i;
	    xip = xi + le;
	    temp.real = xi->real + xip->real;
	    temp.imag = xi->imag + xip->imag;
	    xip->real = xi->real - xip->real;
	    xip->imag = xi->imag - xip->imag;
	    *xi = temp;
	}

/* remaining iterations use stored w */

	wptr = w + windex - 1;
	for (j = 1 ; j < le ; j++) {
	    u = *wptr;
	    for (i = j ; i < n ; i = i + 2*le) {
		xi = x + i;
		xip = xi + le;
		temp.real = xi->real + xip->real;
		temp.imag = xi->imag + xip->imag;
		tm.real = xi->real - xip->real;
		tm.imag = xi->imag - xip->imag;
		xip->real = tm.real*u.real - tm.imag*u.imag;
		xip->imag = tm.real*u.imag + tm.imag*u.real;
		*xi = temp;
	    }
	    wptr = wptr + windex;
	}
	windex = 2*windex;
    }

/* rearrange data by bit reversing */

    j = 0;
    for (i = 1 ; i < (n-1) ; i++) {
	k = n/2;
	while(k <= j) {
	    j = j - k;
	    k = k/2;
	}
	j = j + k;
	if (i < j) {
	    xi = x + i;
	    xj = x + j;
	    temp = *xj;
	    *xj = *xi;
	    *xi = temp;
	}
    }
}

/**************************************************************************
log2 - 基2求对数函数

函数返回以2为底的输入整数的对数值的整数.
如果对数值位于两个整数之间,则返回较大值.

int log2(unsigned int x)
*************************************************************************/

int log2(x)
    unsigned int x;
{
    unsigned int mask,i;

    if(x == 0) return(-1);      /* zero is an error, return -1 */

    x--;                        /* get the max index, x-1 */

    for(mask = 1 , i = 0 ; ; mask *= 2 , i++) {
	if(x == 0) return(i);   /* return log2 if all zero */
	x = x & (~mask);        /* AND off a bit */
    }
}

/**************************************************************************
ADDNOISE  - 在输入信号上加入均匀分布或高斯分布的噪声

输入参数:
	 double  *Indata: 输入的未加噪声的信号序列;
	 unsigned length: 输入信号序列的长度;
	 int      type  : 加噪声的类型, 0:均匀  1:高斯
	 double   nmult : 加入噪声的幅度值系数;

输出参数:
	 double  *Indata: 在输入序列上加入噪声后,仍然将数据
			  放在 *Indata 数据区中.
*************************************************************************/

double *addnoise(double *Indata,unsigned length,int type,double nmult)
{
    int i,j;
    double gaussian(),uniform();

/* add noise to record */
    if(type == 0)
	for(i = 0 ; i < length ; i++) Indata[i] += nmult*uniform();
    else
	for(i = 0 ; i < length ; i++) Indata[i] += nmult*gaussian();
    return(Indata);
}

/**************************************************************************
gaussian - 产生均值为零,单位方差的高斯分布随机数

函数返回 double型零均值,单位方差的高斯分布随机数.
使用 Box-Muller 转换方法,将一对均匀分布的随机变量映射为一对
高斯随机变量.

double gaussian()
*************************************************************************/

double gaussian(void)
{