IIR滤波器的设计
1. P1~P7为.h文件中的程序
2. P8为main.c文件中的程序
3. P9为MATLAB程序
4. 最后是解释的部分
5. tmwtypes.h不是本次作业的重点，不放在本文档中
MATLAB fdatool生成的.h文件如下：
/*
* Filter Coefficients (C Source) generated by the Filter Design and Analysis Tool
* Generated by MATLAB(R) 9.3 and Signal Processing Toolbox 7.5.
* Generated on: 01-May-2018 19:51:08
*/
/*
* Discrete-Time IIR Filter (real)
* -------------------------------
* Filter Structure : Direct-Form II, Second-Order Sections
* Number of Sections : 19
* Stable : Yes
* Linear Phase : No
*/
/* General type conversion for MATLAB generated C-code */
#include "tmwtypes.h"
/*
* Expected path to tmwtypes.h
* D:\Matlab2017b\extern\include\tmwtypes.h
*/
#define MWSPT_NSEC 39
const int B_LEN[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,2,1 }; //这个是下面B数组的每一行包含的实际有效的元素的个数
real64_T B[MWSPT_NSEC][3] = { //这个就是IIR滤波器传递函数里的B
{
0.147266225098, 0, 0 //这个是第一级的放大倍数（第一级分子的放大倍数）
},
{
1, 2, 1 //这个是第一级的b0,b1,b2,下同
},
{
0.1390753664169, 0, 0
},
{
1, 2, 1
},
{
0.1318227971961, 0, 0
},
{
1, 2, 1
},
{
0.1254019296567, 0, 0
},
{
1, 2, 1
},
{
0.1197203276969, 0, 0
},
{
1, 2, 1
},
{
0.1146981171376, 0, 0
},
{
1, 2, 1
},
{
0.1102664769354, 0, 0
},
{
1, 2, 1
},
{
0.1063662631503, 0, 0
},
{
1, 2, 1
},
{
0.1029467861143, 0, 0
},
{
1, 2, 1
},
{
0.09996474328992, 0, 0
},
{
1, 2, 1
},
{
0.09738330065735, 0, 0
},
{
1, 2, 1
},
{
0.09517131086323, 0, 0
},
{
1, 2, 1
},
{
0.09330265475678, 0, 0
},
{
1, 2, 1
},
{
0.09175569304301, 0, 0
},
{
1, 2, 1
},
{
0.09051281581858, 0, 0
},
{
1, 2, 1
},
{
0.08956007925443, 0, 0
},
{
1, 2, 1
},
{
0.08888692038072, 0, 0
},
{
1, 2, 1
},
{
0.08848594266561, 0, 0
},
{
1, 2, 1
},
{
0.2972419330929, 0, 0
},
{
1, 1, 0
},
{
1, 0, 0 //这个是滤波器输出分子总的放大倍数
}
};
const int A_LEN[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,
3,1,2,1 }; //这个是下面A数组的每一行包含的实际有效的元素的个数
real64_T A[MWSPT_NSEC][3] = {
{
1, 0, 0 //第一级的放大倍数（第一级的分母的放大倍数）
},
{
1, -1.351827054515, 0.9408919549076 //第一级的a0,a1,a2，a0总是1，下同
},
{
1, 0, 0
},
{
1, -1.276639248503, 0.8329407141705
},
{
1, 0, 0
},
{
1, -1.210064449829, 0.7373556386132
},
{
1, 0, 0
},
{
1, -1.151124238334, 0.6527319569608
},
{
1, 0, 0
},
{
1, -1.09897009887, 0.5778514096574
},
{
1, 0, 0
},
{
1, -1.052868828174, 0.5116612967244
},
{
1, 0, 0
},
{
1, -1.012188685003, 0.4532545927443
},
{
1, 0, 0
},
{
1, -0.9763867588688, 0.40185181147
},
{
1, 0, 0
},
{
1, -0.9449977450844, 0.3567848895416
},
{
1, 0, 0
},
{
1, -0.9176241489659, 0.3174831221256
},
{
1, 0, 0
},
{
1, -0.8939278534435, 0.2834610560729
},
{
1, 0, 0
},
{
1, -0.8736229420763, 0.2543081855292
},
{
1, 0, 0
},
{
1, -0.8564696546975, 0.2296802737247
},
{
1, 0, 0
},
{
1, -0.8422693538779, 0.2092921260499
},
{
1, 0, 0
},
{
1, -0.8308603898991, 0.1929116531734
},
{
1, 0, 0
},
{
1, -0.8221147656912, 0.1803550827089
},
{
1, 0, 0
},
{
1, -0.8159355187059, 0.1714832002288
},
{
1, 0, 0
},
{
1, -0.8122547526431, 0.1661985233055
},
{
1, 0, 0
},
{
1, -0.4055161338143, 0
},
{
1, 0, 0 //这个是滤波器输出分母总的放大倍数
}
};
Main.c文件如下：
#include<stdio.h>
#include<math.h>
#include"iir0.h"
#define IIR 19 //节数
#define pi 3.1415926
#define fs 48000 //采样频率 fs
#define N 128 //时域采样点数
#define f0 3000 //频率1
#define f1 12000 //频率2
void iir()
{
FILE *fp;
int i;
double x[N];
fp=fopen("E://I love DSP//work2.txt", "w");
double w[19][3]={};
for (i=0;i<N;i++)
{
int k;
double temp;
x[i]=cos(2*pi*f0/fs*i)+cos(2*pi*f1/fs*i); //被滤信号
temp=x[i];
for(k=0;k<IIR;k++)
{
w[k][0]=temp-A[2*k+1][1]*w[k][1]-A[2*k+1][2]*w[k][2];
temp=B[2*k+1][0]*w[k][0]+B[2*k+1][1]*w[k][1]+B[2*k+1][2]*w[k][2];
w[k][2]=w[k][1];
w[k][1]=w[k][0];
temp*= B[2*k][0]; //增益不为1，不可忽略
}
fprintf(fp, "%lf\n", temp);
}
fclose(fp);
}
int main()
{
iir();
return 0;
}
Matlab程序如下：
clf;
fp=fopen("E://I love DSP//work2.txt","r");
x = fscanf(fp,'%f');
fs=48000;N=128; %采样频率和数据点数
n=0:N-1;
t=n/fs; %时间序列
y=fft(x,N); %对信号进行快速Fourier变换
mag=abs(y); %求得Fourier变换后的振幅
f=n*fs/N; %频率序列
stem(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=100');grid on;
MATLAB生成的fft图像为：
由C程序可见,待滤波的数据有3000Hz和12000Hz两种频率
由MATLAB绘制的频谱图可见，除了少许谐波外，主要振幅最大的频率为3000Hz，12000Hz已经不存在，滤波成功
9165
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
