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已經不存在，濾波成功