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文件程序如下所示：

clc
N=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滤波并FFT
x = 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结果：

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

记录用matlab的FDA工具箱产生系数的使用，免得每次都要想一想

maltlab的filter文档里介绍了滤波器实现方式如下

 example:

             用FDA产生二阶高通IIR,参数如下

，系数为a = [1 -1.1429 0.4128]，b = [1 2 -1],用matlab代码实现如下；

t=0:0.0001:1;  
s=sin(2*pi*100*t)+sin(2*pi*2500*t)+cos(2*pi*800*t);
figure,plot(abs(fft(s))),title('fft');
b = [1 -2 1];
a = [1 -1.1429 0.4128];
s2 = filter(b,a,s);
s3 = s;
s3(1) = b(1)*s(1);
s3(2) = -1*a(2)*s(1)+b(1)*s(2)+b(2)*s(1);
for i = 3:10000
    s3(i) = -1*(a(1+1)*s3(i-1)+a(1+2)*s3(i-2))+b(0+1)*s(i)+b(1+1)*s(i-1)+b(2+1)*s(i-2);
end
figure,plot(abs(fft(s2))) ,title('filter');
figure,plot(abs(fft(s3))),title('diff')

容易弄错的地方，matlab是以以下形式给出的系数，对应的z变化形式分母为

对应的差分方程为