C数字滤波器

复习数字信号处理，根据IIR滤波器的差分方程写了段相应的C程序，滤波器系数用matlab产生，滤波程序如下：

void IIR(const double *x,double *y){    //滤波器系数,低通，截止频率1KHz    double B[3]={0,0.116,0.0715};     //分子系数，系统零点    double A[4]={1,1.6 ,-1.02 ,0.232};//分母系数，系统极点int n=0;y[0]=B[0]*x[0];y[1]=A[1]*y[0]+B[0]*x[1]+B[1]*x[0];y[2]=A[1]*y[1]+A[2]*y[0]+B[0]*x[2]+B[1]*x[1]+B[2]*x[0];for(n=3;n<N;n++)y[n]=A[1]*y[n-1]+A[2]*y[n-2]+A[3]*y[n-3]+B[0]*x[n]+B[1]*x[n-1]+B[2]*x[n-2];}

用matlab产生信号并写入到文件中，在vs2012里用c/c++读取文件数据进行滤波并做FFT，并将结果写入到文件，最后用matlab读取结果文件中的数据绘图，与直接用matlab滤波及FFT做对比，其中C与matlab都是对数据做2048点FFT,C的FFT程序之前已调试成功，这里只将其移植调用。matlab文件程序如下所示：

clcN=2048;%向文件中写入原始数据t=0:0.0001:1;s=sin(2*pi*100*t)+sin(2*pi*2500*t)+cos(2*pi*800*t);             %??0~pi/2?Sin?;fidc=fopen('D:\Microsoft Visual Studio 11.0\Projects\TestC++\DIT-FFT\DIT-FFT\InData.txt','wt');       %?"wt"?????,\n??? for t=1:10000; fprintf(fidc,'%f,', s(t)); end fclose(fidc);  %从文件读入数据并显示频谱x = importdata('D:\Microsoft Visual Studio 11.0\Projects\TestC++\DIT-FFT\DIT-FFT\InData.txt');y=fft(x*10000/N,N);y=abs(y(1:N/2));figure(1)subplot(1,2,1);plot(y)%读入C程序FFT结果并显示x_fft = importdata('D:\Microsoft Visual Studio 11.0\Projects\TestC++\DIT-FFT\DIT-FFT\InData_FFT.txt');subplot(1,2,2);plot(x_fft(1:N/2));%从文件读入数据，经过matlab滤波并FFTx = importdata('D:\Microsoft Visual Studio 11.0\Projects\TestC++\DIT-FFT\DIT-FFT\InData.txt');z=[0 0.116 0.0715];p=[1 -1.6 1.02 -0.232];y2=filter(z,p,x);y3=fft(y2,N);y3=abs(y3(1:N/2));figure(2);subplot(1,2,1);plot(y3);%从文件读入经过C滤波并FFT的结果x_iir_fft = importdata('D:\Microsoft Visual Studio 11.0\Projects\TestC++\DIT-FFT\DIT-FFT\OutData_FFT.txt');subplot(1,2,2);plot(x_iir_fft(1:N/2));

matlab的FFT与c程序的FFT:

MATLAB滤波及FFT结果与C滤波及FFT结果：

结果比较程序大致运行正确。

程序目前目前比较粗糙，还有很多地方有待优化，算法原理理解还不是很透彻，待日后慢慢研究……