C++實作matlab的fir1函數

• 導言
• • 函數需求分析
  • 數學過程
  • 源碼
  • 注意

導言

最近在進行Qt開發，涉及大量的matlab轉C的工作，其中包括插值濾波等，但遺憾就求濾波系數的函數fir1而言，多數都是直接用的matlab生成的系數進行濾波，很少有用C生成的，有少數用C進行實作也和matlab生成的系數相差甚遠，是以這裡對matlab生成濾波系數的fir1函數進行了研究并用C++進行了重寫，以備進行FIR濾波。

經對比，結果與matlab的fir1函數生成的濾波系數完全一緻。

函數需求分析

這裡我對工作中調用matlab的fir1()函數的執行個體為例子進行講解：

1. 原matlab程式 ：由256（fir1()函數裡進行了+1）階，截止頻率為0.25的低通濾波器（預設使用hamming窗）

lbn=255;
   lbf=fir1(lbn,0.25,'low');
           

1. fir1（）生成的結果

   

2. 分析

   可以看到，輸入參數為階數n，截止頻率Wn，選擇低通濾波（也是預設的選擇），輸出參數為一個一維數組，即濾波參數。

3. C++版fir1（）輸入輸出參數詳情：

參數名	介紹
n	階數，對應matlab的fir1裡的階數n
Wn	對應matlab的fir1裡的階數Wn，但應注意傳進來的資料應存在一個vector的double數組裡。
h	double類型的vector數組，裡面存的是濾波器系數

vector <double> fir1(int numtaps, vector<double> cutoff)
  /*
  	未寫排錯  檢查輸入有需要自己進行完善
  	原matlab函數fir(n, wn)	【函數預設使用hamming】

  	參數輸入介紹：
  		n：  對應matlab的fir1裡的階數n
  		Wn:  對應matlab的fir1裡的階數Wn，但應注意傳進
  				 來的資料應存在一個vector的double數組裡。

  	參數輸出介紹：
  				vector <double>的一個數組，裡面存的是長度
  				為n的濾波器系數。
  */;
           

1. C++版fir1函數生成結果 ：結果表明，生成結果與matlab生成結果一緻。

   

數學過程

數學小辣雞的我，希望得到大家的指正

源碼

那麼廢話不多說，直接上代碼

#include <stdio.h>
#include <math.h>
#include <vector>
#include <iostream>
using namespace std;
#define PI acos(-1)
vector<double> sinc(vector<double> x)
{
   vector<double> y;
   for (int i = 0; i < x.size(); i++)
   {
   	double temp = PI * x[i];
   	if (temp == 0) {
   		y.push_back(0.0);
   	}
   	else {
   		y.push_back(sin(temp) / temp);
   	}
   }
   return y;
}

vector <double> fir1(int n, vector<double> Wn)
{

   /*
   	未寫排錯  檢查輸入有需要自己進行完善
   	原matlab函數fir(n, wn)	【函數預設使用hamming】

   	參數輸入介紹：
   		n：  對應matlab的fir1裡的階數n
   		Wn:  對應matlab的fir1裡的階數Wn，但應注意傳進
   				 來的資料應存在一個vector的double數組裡。

   	參數輸出介紹：
   				vector <double>的一個數組，裡面存的是長度
   				為n的濾波器系數。
   */
   
   //在截止點的左端插入0（在右邊插入1也行）
   //使得截止點的長度為偶數，并且每一對截止點對應于通帶。
   if (Wn.size() == 1 || Wn.size() % 2 != 0) {
   	Wn.insert(Wn.begin(), 0.0);
   }

   /*
   	‘ bands’是一個二維數組，每行給出一個 passband 的左右邊緣。
   	（即每2個元素組成一個區間）
   */
   vector<vector <double>> bands;
   for (int i = 0; i < Wn.size();) {
   	vector<double> temp = { Wn[i], Wn[i + 1] };
   	bands.push_back(temp);
   	i = i + 2;
   }

   // 建立系數
   /*
   	m = [0-(n-1)/2,
   		 1-(n-1)/2,
   		 2-(n-1)/2,
   		 ......
   		 255-(n-1)/2]
   	h = [0,0,0......,0]
   */
   double alpha = 0.5 * (n - 1);
   vector<double> m;
   vector<double> h;
   for (int i = 0; i < n; i++) {
   	m.push_back(i - alpha);
   	h.push_back(0);
   }
   /*
   	對于一組區間的h計算
   	left:	一組區間的左邊界
   	right:  一組區間的右邊界
   */
   for (int i = 0; i < Wn.size();) {
   	double left = Wn[i];
   	double right = Wn[i+1];
   	vector<double> R_sin, L_sin;
   	for (int j = 0; j < m.size(); j++) {
   		R_sin.push_back(right * m[j]);
   		L_sin.push_back(left * m[j]);
   	}
   	for (int j = 0; j < R_sin.size(); j++) {
   		h[j] += right * sinc(R_sin)[j];
   		h[j] -= left * sinc(L_sin)[j];
   	}

   	i = i + 2;
   }

   // 應用視窗函數，這裡和matlab一樣
   // 預設使用hamming，要用别的窗可以去matlab查對應窗的公式。
   vector <double> Win;
   for (int i = 0; i < n; i++)
   {
   	Win.push_back(0.54 - 0.46*cos(2.0 * PI * i / (n - 1)));	//hamming窗系數計算公式
   	h[i] *= Win[i];
   }

   bool scale = TRUE;
   // 如果需要，現在可以處理縮放.
   if (scale) {
   	double left = bands[0][0];
   	double right = bands[0][1];
   	double scale_frequency = 0.0;
   	if (left == 0)
   		scale_frequency = 0.0;
   	else if (right == 1)
   		scale_frequency = 1.0;
   	else
   		scale_frequency = 0.5 * (left + right);

   	vector<double> c;
   	for (int i = 0; i < m.size(); i++) {
   		c.push_back(cos(PI * m[i] * scale_frequency));
   	}
   	double s = 0.0;
   	for (int i = 0; i < h.size(); i++) {
   		s += h[i] * c[i];
   	}
   	for (int i = 0; i < h.size(); i++) {
   		h[i] /= s;
   	}
   }
   return h;
}
           

注意

該文章僅個人學習使用，歡迎大家一起交流學習