频率采样法设计fir滤波器程序 python

频率采样法FIR滤波器设计

• 频率采样法
• code

上一篇中，简单介绍了FIR滤波器的分类以及窗函数设计低通滤波器的方法，这里总结一下另一常用的频率采样法

频率采样法

窗函数法是从时域的角度出发，把理想的非因果无限长的单位脉冲响应$h_d(n)$截断为因果有限长的$h(n)$，
　　而频率采样法，是直接从频率出发，假设咱们有一个目标的频率响应$H_d(k)=H_d(e^{j\omega})|_{\omega_k=2\pi k/N} =H_d(e^{j2\pi k/N})$
　　这个公式的意思是，咱们有个目标频率响应$H_d(e^{j\omega})$，这还是一个关于$\omega$的连续函数，现在在$\left | e^{j\omega}\right |=1$的单位圆上等间隔采样N个点，得到$H_d(k)$，这个就是咱们实际得到的频率响应，最后，傅里叶反变换到时域，就得到了咱们需要的脉冲响应。
　　$h(n) = \frac{1}{N}\sum_{k=0}^{N-1}H_d(k)e^{j2\pi kn/N},n=0,1,....N-1$
　　当然，一般在时域给出需要的频率响应$H_d(k)$的时候，一般只会给出$0 \sim \pi$的值，那么在做傅里叶反变换前就需要将$H_d(k)$的$\pi \sim 2\pi$的值按FIR滤波器线性相位的条件补起来。
　　定义理想的频率响应为
　　$H(e^{j\omega})=H(\omega)e^{j\theta(\omega)}$
　　表示为幅度响应和相位响应相乘
　　这里的满足线性相位条件可以表示为如下：
　　

其中，
$\theta(k)=-\frac{N-1}{N}\pi k$
上面的中括号表示取整，当N为偶数时，$H(\frac{N}{2})=0$
其实这个条件就是需要频率响应满足共轭对称的性质，跟DFT中的一样，可以写成下面这样
$H_{d}(N-k)=H_{d}^{*}(k),k=1,2,....,[\frac{N-1}{2}]$

code

测试下冲击响应分别为奇偶的情况

%   线性相位FIR,h(n)为实序列，且h(n) = ±h(N-1-n),即h(n)对称（奇对称或偶对称）
%   N为奇数，h(n)偶对称，幅度A（w)满足偶对称
%                           A(k) = A(N-k) ,k = 0,1...N-1
%                       相位phy = -w*(N-1)/2
%                                                 = -2*pi*k/N*(N-1)/2
%                                                 = -k*pi*(N-1)/N
%   N为偶数，h(n)偶对称，幅度A(w)满足奇对称
%                           A(k) = -A(N-k) ,k = 0,1...N-1
%                       相位phy = -k*pi*(N-1)/N
%
%   H(w) = A(w)*exp(j*phy)
%        = A(w)*exp(-j*k*pi*(N-1)/N)  最终H(w)满足共轭对称性
%
%
%
%%
close all
N = 33;    % N为奇数
fs = 16000;
fc = 4000.0;               %cut-off frequency
wc = fix(fc*N/fs)+1;


Hd = zeros(1,(N+1)/2);
Hd(1) = 1;  % DC
Hd(2:wc-1) = 1; %passband
Hd(wc) = 0.4;% transition attenuation
Hd(wc+1:end) = 0;   % stopband attenuation

k = 0:(N-1)/2;
A = exp(-1j*pi*k*(N-1)/N);

Hd = Hd.*A;

Hk = [Hd,conj(fliplr(Hd(2:end)))];


h1 = ifft(Hk);
figure,freqz(h1)

%%
N = 32;    % N为偶数
fs = 16000;
fc = 4000.0;               %cut-off frequency
wc = fix(fc*N/fs)+1;


Hd = zeros(1,N/2+1);
Hd(1) = 1;  % DC
Hd(2:wc-1) = 1; %passband
Hd(wc) = 0.4;% transition attenuation
Hd(wc+1:end-1) = 0;   % stopband attenuation
Hd(N/2) = 0;  % Nyquist frequency


k = 0:N/2;
A = exp(-1j*pi*k*(N-1)/N);

Hd = Hd.*A;

Hk = [Hd,conj(fliplr(Hd(2:end-1)))];


h2 = ifft(Hk);
figure,freqz(h2)

可以看到得到的频率响应以及时域的冲击响应如下