使用C++來寫一個IIR濾波器
我們首先要在MATLAB中設計一個IIR濾波器，並生成一個頭文件，這個頭文件中反映了IIR濾波器的頻率響應特性
理論支持
IIR濾波叫做遞歸濾波器，它是一種具有反饋的濾波器。當階數較大時一般采取多個二階節濾波進行串聯，這樣可以提高系統穩定性。
一個二階節系數規律如圖所示：
可以寫出第K個二階節的差分方程
N個二階節的級聯結構如下圖所示：
根據二階節圖，把前一級的輸出作為后一級的輸入，就可以通過軟件實現IIR數字濾波的功能。
使用Matlab生成頭文件
首先打開MATLAB中Filter Design & Analysis Tool
這里我們先設計一個低通濾波器
Fs代表采樣頻率，采樣頻率必須大於原信號最高頻率的兩倍，
否則會產生頻譜混疊。
Fpass為通帶頻率，Fstop為阻帶截止頻率
這些參數設置好就可以點擊Design Filter
生成的是一個二階節濾波組合，一共有31階，也就是多個二階濾波器的組合
接着在Target選項中生成C Header File
Numerator為分子系數數組的命名，Numerator length為分子系數數組的長度，
Denominator為分母。
對生成頭文件進行分析
以下以Fpass為10K，Fstop為12K的低通濾波器舉栵
在使用頭文件前需要根據情況將Matlab的數據類型轉換為C++支持的數據類型，這里我們使用double類型
在分析頭文件前先看下Matlab提供的第一節濾波參數
以第一個二階節的數據舉例：
- Numerator: 1 2 1
- Denominator: 1 -0.55930961405944157 0.92579835996619642
- Gain：0.34162218647668868
Numerator為分子的系數，分別為b0，b1，b2
Denominator為分母的系數，分別為a0，a1，a2.
Gain為各節的增益，此項為了穩定各節，穩定信號大小
接着對照頭文件，以下為頭文件主要部分的一段截取：
#define MWSPT_NSEC 41 int NL[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,3,1,2,1 }; double NUM[MWSPT_NSEC][3] = { { 0.3416221864767, 0, 0 }, { 1, 2, 1 }, { 0.3180955154747, 0, 0 }, { 1, 2, 1 },......
MWSPT_NSEC為濾波器階數，具體的節數在頭文件開頭的注釋中
NL[MWSPT_NSEC]這個數組定義了NUM[MWSPT_NSEC][3]數組每一行的有用數據個數（可以不用）
在NUM[MWSPT_NSEC][3]數組（分子參數）奇數行第一項都為增益項，偶數行為3個系數，分別為b0,b1,b2。
由此可以找出規律，定義K為目前所在的階數，p為數組的首指針，則，每一節的增益項為（p+6*K）,第一個系數為（p+3+6*K），
第二個系數為（p+3+6*K+1），第三個系數為（p+3+6*K+2）。
C++編程實現
在軟件設計的過程中，每個二階節的延遲變量只取 和 ， 作為中間變量在過程中直接賦給 。這是因為對於下一個輸入數據n+1的延遲變量即為上一個輸入數據的 和 ，采用這種方式進行設計，可以節省寄存器的空間。
為了提高處理速度，程序中需要使用指針進行參數傳遞，特別注意二維數組的首地址傳遞方式為&a[0][0]->double* a。
濾波函數
double iir(double *a, double *b,double* w, double xin, int N_IIR)
{
int k;
double temp = xin;
for (k = 0; k<N_IIR; k++)
{
*(w+k*3) = temp - *(a + 3+6 * k + 1) *(*(w + k * 3+1)) - *(a + 3 + 6 * k + 2) *(*(w + k * 3+2));
//這里temp為本二階節的輸入，也是上一個二階節的輸出
temp = *(b + 3 + 6 * k )* (*(w + k * 3)) + *(b + 3 + 6 * k + 1) * (*(w + k * 3+1)) + *(b + 3+6 * k + 2)* (*(w + k * 3+2));
//這里temp為本二階節的輸出，也是下一個二階節的輸入
*(w + k * 3 + 2) = *(w + k * 3 + 1);
*(w + k * 3 + 1) = *(w + k * 3);
temp = temp*(*(b + 6 * k));//放大倍數，穩定信號
}
return temp;
}
實際測試
測試Fpass為10K，Fstop為12K的低通濾波器
在程序中輸入三個頻率為2K，11K，20K的信號，理論上2k完全通過，11k部分衰減，20K完全濾除。
上圖為原信號，下圖為濾波后信號。
實際測試發現符合設計要求，而且在過渡帶信號也基本完全衰減。
測試用C++程序 void main() { const int N = 100; int i,j; double xn[N]; double w[20][3]; double yn[N]; for (i = 0; i < 20; i++)//初始化 { for (j = 0; j < 3; j++) w[i][j] = 0; } for (i = 0; i < N; i++) { xn[i] = sin(2 * 3.1416 * 20 / 50 * i)+ sin(2 * 3.1416 * 2 / 50 * i)+ sin(2 * 3.1416 * 11 / 50 * i); yn[i]=iir(&DEN[0][0], &NUM[0][0], &w[0][0],xn[i], 20); } ofstream SaveFile_a("xn.txt"); for (i = 0; i<N; i++) SaveFile_a << " " << xn[i] << endl; SaveFile_a.close(); ofstream SaveFile_b("yn.txt"); for (i = 0; i<N; i++) SaveFile_b << " " << yn[i] << endl; SaveFile_a.close(); } 分析用Matlab程序 xn1=fopen('xn.txt','r'); [xn,count]=fscanf(xn1,'%f'); fclose(xn1); N = length(xn);%求取抽樣點數 xn_f = fft(xn);%對信號進行傅里葉變換 xn_f=abs(xn_f(1:N/2)); f = 50000/N*(0:N/2-1); subplot(211); stem(f,abs(xn_f)); xlabel('Frequency / (s)');ylabel('Amplitude'); title('原信號頻譜'); grid; yn1=fopen('yn.txt','r'); [yn,count]=fscanf(yn1,'%f'); fclose(yn1); yn_f = fft(yn);%對信號進行傅里葉變換 yn_f=abs(yn_f(1:N/2)); subplot(212); stem(f,abs(yn_f)); xlabel('Frequency / (s)');ylabel('Amplitude'); title('濾波后信號頻譜'); grid;