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

其中,
θ
(
k
)
=
−
N
−
1
N
π
k
\theta(k)=-\frac{N-1}{N}\pi k
θ(k)=−NN−1πk
上面的中括号表示取整,当N为偶数时,
H
(
N
2
)
=
0
H(\frac{N}{2})=0
H(2N)=0
其实这个条件就是需要频率响应满足共轭对称的性质,跟DFT中的一样,可以写成下面这样
H
d
(
N
−
k
)
=
H
d
∗
(
k
)
,
k
=
1
,
2
,
.
.
.
.
,
[
N
−
1
2
]
H_{d}(N-k)=H_{d}^{*}(k),k=1,2,....,[\frac{N-1}{2}]
Hd(N−k)=Hd∗(k),k=1,2,....,[2N−1]
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)
可以看到得到的频率响应以及时域的冲击响应如下