题目
如图1所示,由
x
(
t
)
x(t)
x(t)信号
x
1
(
t
)
x_1(t)
x1(t)、信号
x
2
(
t
)
x_2(t)
x2(t)和信号
x
3
(
t
)
x_3(t)
x3(t)的和构成。对
x
(
t
)
x(t)
x(t)进行均匀采样(采样周期为
T
s
=
0.25
m
s
T_s=0.25ms
Ts=0.25ms),获得其样本信号
x
(
n
T
s
)
x(nT_s)
x(nTs)。
其中,
x
i
(
t
)
=
A
i
c
o
s
(
2
π
f
i
t
)
x_i(t)=A_icos(2\pi f_it)
xi(t)=Aicos(2πfit),参数
A
i
A_i
Ai和
f
i
f_i
fi如下表所示
i= 1 | i= 2 | i= 3 | |
---|---|---|---|
A i A_i Ai | 0.2 | 2 | 0.4 |
f i f_i fi | 0.99kHz | 0.98kHz | 1.6kHz |
问题
(供参考)
傅里叶变换
(1)求 x ( t ) x(t) x(t)的傅立叶变换,并画出其幅频响应曲线。
解:由于 x ( t ) x(t) x(t)由信号 x 1 ( t ) x_1(t) x1(t)、信号 x 2 ( t ) x_2(t) x2(t)和信号 x 3 ( t ) x_3(t) x3(t)的和构成。由此可将 x ( t ) x(t) x(t)表示为:
x ( t ) = ∑ i = 1 3 x i ( t ) = ∑ i = 1 3 A i cos ( 2 π f i t ) x(t)= \sum\limits_{i = 1}^3 {{x_i}\left( t \right)} {\rm{ = }}\sum\limits_{i = 1}^3 {{A_i}\cos ( {2\pi {f_i}t} )} x(t)=i=1∑3xi(t)=i=1∑3Aicos(2πfit)
由傅里叶变换公式 X ( j ω ) = ∫ − ∞ + ∞ x ( t ) e − j ω t d t X(j\omega ) = \int_{ - \infty }^{ + \infty } {x(t)} {e^{ - j\omega t}}dt X(jω)=∫−∞+∞x(t)e−jωtdt
F [ x ( t ) ] = F [ ∑ i = 1 3 x i ( t ) ] = ∑ i = 1 3 A i π [ δ ( ω − 2 π f i ) + δ ( ω + 2 π f i ) ] = 0.2 π [ δ ( ω − 2 ∗ 990 π ) + δ ( ω + 2 ∗ 990 π ) ] + 2 π [ δ ( ω − 2 ∗ 980 π ) + δ ( ω + 2 ∗ 980 π ) ] + 0.4 π [ δ ( ω − 2 ∗ 1600 π ) + δ ( ω + 2 ∗ 1600 π ) ] {\cal F}\left[ {x\left( t \right)} \right] = {\cal F}\left[ {\sum\limits_{i = 1}^3 {{x_i}\left( t \right)} } \right]\ = \sum\limits_{i = 1}^3 {{A_i}\pi \left[ {\delta (\omega - 2\pi {f_i}) + \delta (\omega + 2\pi {f_i})} \right]} \\ = 0.2\pi \left[ {\delta (\omega - 2*990\pi ) + \delta (\omega + 2*990\pi )} \right] + \\ 2\pi \left[ {\delta (\omega - 2*980\pi ) + \delta (\omega + 2*980\pi )} \right] + \\ 0.4\pi \left[ {\delta (\omega - 2*1600\pi ) + \delta (\omega + 2*1600\pi )} \right] F[x(t)]=F[i=1∑3xi(t)] =i=1∑3Aiπ[δ(ω−2πfi)+δ(ω+2πfi)]=0.2π[δ(ω−2∗990π)+δ(ω+2∗990π)]+2π[δ(ω−2∗980π)+δ(ω+2∗980π)]+0.4π[δ(ω−2∗1600π)+δ(ω+2∗1600π)]
%傅里叶变换
clc;
clear;
syms t;
x1=0.2*cos(2*pi*1280*t);
x2=2*cos(2*pi*1290*t);
x3=0.4*cos(2*pi*1800*t);
x=x1+x2+x3; %x(t)信号由x1、x2、x3的和组成。
F=fourier(x)
%幅频响应曲线
close all;
clear;
clc;
Ts=0.00025; %采样周期,以0.25ms的周期均匀采样
Fs=1/Ts; %采样频率
t=0:Ts:1-Ts; %采样间隔时长
x1=0.2*cos(2*pi*990*t); %三分量信号
x2=2*cos(2*pi*980*t);
x3=0.4*cos(2*pi*1600*t);
x=x1+x2+x3; %x(t)信号由x1、x2、x3的和组成。
figure(1);
subplot(2,1,1);
plot(t,x);
xlabel('时间/s');ylabel('幅值');
title("original signal");%绘制原始信号图
%频域
N=(1-Ts)/Ts+1;%采样点数
n=0:N-1;
fft_x=n*Fs/N-Fs/2;
x_fft=fft(x);
shift_f=abs(fftshift(x_fft));
y2=shift_f/N;
subplot(2,1,2);
plot(fft_x,y2);
xlabel('频率/Hz');ylabel('幅值');
title('FFT'); %0频率分量位于坐标中心
grid on;
x
3
(
t
)
x_3(t)
x3(t)信号:
频谱图:
最少采样点数
(2)如果希望能从所获取到的样本信号
x
(
n
T
s
)
x(nT_s)
x(nTs)中,通过DFT有效的区分出信号
x
1
(
t
)
x_1(t)
x1(t)、
x
2
(
t
)
x_2(t)
x2(t)和
x
3
(
t
)
x_3(t)
x3(t),试求最少需要的采样点数N,画出此时样本信号
x
(
n
T
s
)
x(nT_s)
x(nTs)(2)的离散傅里叶变换(DFT)的主值区间的频谱图,并标出数字频率及其对应的模拟频率。
解:分析可知原时域模拟信号为
x
(
t
)
=
∑
i
=
1
3
x
i
(
t
)
=
∑
i
=
1
3
A
i
cos
(
2
π
f
i
t
)
=
0.2
cos
(
2
π
×
990
t
)
+
2
cos
(
2
π
×
980
t
)
+
0.4
cos
(
2
π
×
1600
t
)
x\left( t \right) = \sum\limits_{i = 1}^3 {{x_i}\left( t \right)} {\rm{ = }}\sum\limits_{i = 1}^3 {{A_i}\cos \left( {2\pi {f_i}t} \right)} = 0.2\cos \left( {2\pi \times 990t} \right) + 2\cos \left( {2\pi \times 980t} \right) + 0.4\cos \left( {2\pi \times 1600t} \right)
x(t)=i=1∑3xi(t)=i=1∑3Aicos(2πfit)=0.2cos(2π×990t)+2cos(2π×980t)+0.4cos(2π×1600t)
即信号的3个频率分量为
f
1
=
990
H
z
f_1=990Hz
f1=990Hz,
f
2
=
980
H
z
f_2=980Hz
f2=980Hz,
f
3
=
1600
H
z
f_3=1600Hz
f3=1600Hz
由采样周期为
T
s
=
0.25
m
s
T_s=0.25ms
Ts=0.25ms,根据抽样间隔与抽样频率的关系可知
f
s
=
1
T
s
=
1
0.25
m
s
=
4000
H
z
{f_s} = \frac{1}{{{T_s}}} = \frac{1}{{0.25ms}} = 4000Hz
fs=Ts1=0.25ms1=4000Hz
抽样后得到的离散序列为
x
(
n
)
=
x
a
(
t
)
∣
t
=
n
T
s
=
0.2
cos
(
99
π
200
n
)
+
2
cos
(
49
π
100
n
)
+
0.4
cos
(
4
π
5
n
)
x(n) = {x_a}(t){|_{t = n{T_s}}} = 0.2\cos (\frac{{99\pi }}{{200}}n) + 2\cos (\frac{{49\pi }}{{100}}n) + 0.4\cos (\frac{{4\pi }}{5}n)
x(n)=xa(t)∣t=nTs=0.2cos(20099πn)+2cos(10049πn)+0.4cos(54πn)
此时三个频率分量
2
π
ω
1
=
2
π
×
200
99
π
=
400
99
2
π
ω
2
=
2
π
×
100
49
π
=
200
49
2
π
ω
1
=
2
π
×
5
4
π
=
5
2
\begin{array}{l} \frac{{2\pi }}{{{\omega _1}}} = 2\pi \times \frac{{200}}{{99\pi }} = \frac{{400}}{{99}}\\ \frac{{2\pi }}{{{\omega _2}}} = 2\pi \times \frac{{100}}{{49\pi }} = \frac{{200}}{{49}}\\ \frac{{2\pi }}{{{\omega _1}}} = 2\pi \times \frac{5}{{4\pi }} = \frac{5}{2} \end{array}
ω12π=2π×99π200=99400ω22π=2π×49π100=49200ω12π=2π×4π5=25
由此可知最小公倍数为N=400,所以x(n)的周期为N=400。
一个周期的数据长度为
T
0
=
N
T
s
=
400
×
0.25
m
s
=
0.1
s
{T_0} = N{T_s} = 400 \times 0.25ms = 0.1s
T0=NTs=400×0.25ms=0.1s
当抽样点数为周期的整数倍时便不会产生频谱泄露。
由于信号中最小的频率间隔为
(
f
1
−
f
2
=
10
H
z
)
({f_1} - {f_2} = 10Hz)
(f1−f2=10Hz)所以
F
0
=
10
H
z
F_0=10Hz
F0=10Hz,为此所需的数据长度为
T
0
=
1
F
0
=
1
10
H
z
=
0.1
s
{T_0} = \frac{1}{{{F_0}}} = \frac{1}{{10Hz}} = 0.1s
T0=F01=10Hz1=0.1s
与一个周期的数据长度相同,可取抽样点数
N
≥
f
s
f
1
−
f
2
=
400
N \ge \frac{{{{\rm{f}}_{\rm{s}}}}}{{{f_1} - {f_2}}} = 400
N≥f1−f2fs=400,因此最少需要的采样点数N取400。
数字角频率
ω
=
2
π
f
f
S
\omega = \frac{{2\pi f}}{{{f_S}}}
ω=fS2πf,得数字频率分别为
ω
1
=
0.154
,
ω
2
=
0.156
,
ω
3
=
0.251
{\omega _1} = 0.154,{\omega _2} = 0.156,{\omega _3} = 0.251
ω1=0.154,ω2=0.156,ω3=0.251。
采用MATLAB方法做频谱分析,结果如图。
clear all;
clc;
Ts=0.00025; %采样间隔为0.25ms
fs=1/Ts;
N=400; %抽样点N取400,抽样频率fs取4000
n=[0:1:N-1];k=[0:1:N-1];
xn=0.2*cos(99*(pi/200)*n)+2*cos(49*(pi/100)*n)+0.4*cos(4*(pi/5)*n);%信号离散序列
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk; %根据定义求信号DFT
X=abs(Xk);
figure(1);
stem(k,X,'*');
grid on;
xlabel('k');
title('幅频特性曲线(K=400,fs=4000Hz)');
ylabel('|X(k)|');
figure(2);
stem(2*pi/fs*k,X,'*');
grid on;
xlabel('数字频率');
title('数字频率(K=400,fs=4000Hz)');
ylabel('|X(k)|');
figure(3);
stem(fs/N*k,X,'*');
grid on;
xlabel('模拟频率');
title('模拟频率(K=400,fs=4000Hz)');
ylabel('|X(k)|');
频谱泄漏
(3)如果采样点数N=512,样本信号 x ( n T s ) x(nT_s) x(nTs)的DFT中,是否存在频谱泄漏,如存在,请分析造成频谱泄露的原因、造成的可能影响,以及可能缓轻它的方法?
解:分析可知当采样点数N为512时,此时的频谱
k
=
125
时
,
f
1
=
k
N
×
f
s
=
125
512
×
4000
=
976.6
H
z
;
k
=
126
时
,
f
1
=
k
N
×
f
s
=
126
512
×
4000
=
984.4
H
z
k = 125时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{125}}{{512}} \times 4000 = 976.6Hz;k = 126时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{126}}{{512}} \times 4000 = 984.4Hz
k=125时,f1=Nk×fs=512125×4000=976.6Hz;k=126时,f1=Nk×fs=512126×4000=984.4Hz此时980Hz的频率分量在此两频率之间。
k
=
126
时,
f
1
=
k
N
×
f
s
=
126
512
×
4000
=
984.4
H
z
;
k
=
127
时,
f
1
=
k
N
×
f
s
=
127
512
×
4000
=
992.2
H
z
k = 126时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{126}}{{512}} \times 4000 = 984.4Hz;k = 127时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{127}}{{512}} \times 4000 = 992.2Hz
k=126时,f1=Nk×fs=512126×4000=984.4Hz;k=127时,f1=Nk×fs=512127×4000=992.2Hz此时990Hz的频率分量在此两频率之间。
k
=
204
时,
f
1
=
k
N
×
f
s
=
204
512
×
4000
=
1593.8
H
z
;
k
=
205
时,
f
1
=
k
N
×
f
s
=
205
512
×
4000
=
1601.6
H
z
k = 204时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{204}}{{512}} \times 4000 = 1593.8Hz;k = 205时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{205}}{{512}} \times 4000 = 1601.6Hz
k=204时,f1=Nk×fs=512204×4000=1593.8Hz;k=205时,f1=Nk×fs=512205×4000=1601.6Hz此时1600Hz的频率分量在此两频率之间。
此时的频率分辨率
F
0
=
f
s
N
=
7.8
<
10
H
z
{F_0} = \frac{{{f_s}}}{N} = 7.8 < 10Hz
F0=Nfs=7.8<10Hz能够满足分辨率的要求,但是会出现少量的频谱泄漏。
造成频谱泄漏的原因:所取的序列点数不是序列周期的整数倍。频谱泄漏的产生与信号时域截断的长度有关,因为在数字信号处理中,信号长度不可能无限长,信号长度越长,分辨率越高,频谱泄漏越少。因为对时域信号进行截断,本质上就是乘了一个矩形窗,也就是相当于在频率域做了卷积。
可能造成的影响:频谱泄露的直接影响是无法准确获得频率幅值,因为在某频率的信号幅值被周围的频率值平分。截短后信号不能准确代表原始信号,谱分析产生误差,出现了不该有的频率分量,从而导致对于频谱的分析不准确。产生的较多的旁瓣还有可能带来谱间干扰和混叠失真。原有谱线被展宽,增加了不必要的频率分量。降低频率分辨率(矩形窗谱主瓣宽度的一半)。谱间干扰,强信号的旁瓣影响弱信号的主瓣(淹没信号)。
缓轻频谱泄漏方法:
1、对周期信号——采样和分析的点数为周期的整倍数;
2、减小主瓣宽度,主瓣宽度与截取长度N成反比,即增加x(n)长度;
3、减小旁瓣宽度,加合适的窗函数,缓慢截断;
4、整周期截断——截取的正好是序列的一个或整数个周期。可以代表原序列,不发生频谱泄漏。
使用MATLAB作采样点N=512时的频谱分析,如图发生了频谱泄漏。
clear all;
clc;
Ts=0.00025; %采样间隔为0.25ms
fs=1/Ts;
N=512; %抽样点N取512,抽样频率fs取4000
n=[0:1:N-1];k=[0:1:N-1];
xn=0.2*cos(99*(pi/200)*n)+2*cos(49*(pi/100)*n)+0.4*cos(4*(pi/5)*n);%信号离散序列
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk; %根据定义求信号DFT
X=abs(Xk);
subplot(2,1,1);
stem(k,X,'*');
grid on;
xlabel('k');
title('幅频特性曲线(K=512,fs=4000Hz)');
ylabel('|X(k)|');
x=fs/N*k;
subplot(2,1,2);
stem(x,X,'*');
grid on;
xlabel('模拟频率');
title('模拟频率(K=512,fs=4000Hz)');
ylabel('|X(k)|');
FIR数字滤波器设计
(4)如希望在样本
x
(
n
T
s
)
x(nT_s)
x(nTs)中,通过FIR数字滤波器有效的抑制
x
3
(
t
)
x_3(t)
x3(t)的频谱分量,使其衰减超过50dB。试设计该数字滤波器(方法不限),给出设计指标及设计方案,并对所设计的滤波器进行仿真,分析其性能。
解:由题意可知,信号
x
1
(
t
)
x_1(t)
x1(t)、
x
2
(
t
)
x_2(t)
x2(t)和
x
3
(
t
)
x_3(t)
x3(t)的频率分别为0.99kHz,0.98kHz和1.6kHz,为有效的抑制
x
3
(
t
)
x_3(t)
x3(t)的频谱分量,使其衰减超过50dB。结合信号频率情况此处令
f
p
=
1.1
k
H
z
,
f
s
t
=
1.5
k
H
z
,
A
s
=
50
d
B
{f_p} = 1.1kHz,{f_{st}} = 1.5kHz,{A_s} = 50dB
fp=1.1kHz,fst=1.5kHz,As=50dB。
考虑采用窗函数设计法。
窗函数法设计线性相位FIR DF的设计步骤为:
1、确定所需滤波器的性能指标要求。
2、按4种理想线性相位DF(低通、带通、带阻,高通)的表达式,将理想滤波器的截止频率(或上下截止频率)确定为过渡带的算术中心频率。
3、得到理想滤波器的单位冲激响应
h
d
(
n
)
h_d(n)
hd(n)。
4、选择窗函数的类型和长度点数N。窗函数的类型由给定滤波器阻带最小衰减As(dB)确定。窗的长度点数N由滤波器过渡带宽来确定。
5、求加窗后的实际滤波器的单位冲激响应
h
(
n
)
=
h
d
(
n
)
ω
(
n
)
h(n) = {h_d}(n)\omega (n)
h(n)=hd(n)ω(n)。
6、检验
H
(
e
j
ω
)
=
D
T
F
T
[
h
(
n
)
]
H({e^{j\omega }}) = DTFT[h(n)]
H(ejω)=DTFT[h(n)]是否满足所求滤波器的性能要求。
将模拟频率转换为数字频率:
ω
p
=
2
π
f
p
/
f
s
=
0.55
π
ω
s
t
=
2
π
f
s
t
/
f
s
=
0.75
π
\begin{array}{l} {\omega _p} = 2\pi {f_p}/{f_s} = 0.55\pi \\ {\omega _{st}} = 2\pi {f_{st}}/{f_s} = 0.75\pi \end{array}
ωp=2πfp/fs=0.55πωst=2πfst/fs=0.75π
理想低通滤波波器的频率响应
H
d
(
e
j
ω
)
{H_d}\left( {{e^{j\omega }}} \right)
Hd(ejω)和单位抽样响应
h
d
(
n
)
{h_d}(n)
hd(n)为:
其中
τ
=
(
N
−
1
)
/
2
\tau = (N - 1)/2
τ=(N−1)/2,
ω
c
=
0.5
ω
p
+
ω
s
t
=
0.65
π
{\omega _c} = 0.5{\omega _p} + {\omega _{st}} = 0.65\pi
ωc=0.5ωp+ωst=0.65π。
由要求的阻带衰减要超过
A
s
=
50
d
B
{A_s} = 50dB
As=50dB,由于海明窗阻带最小衰减为53dB,故选择窗函数的类型为海明窗。
过渡带宽
Δ
ω
=
ω
s
t
−
ω
p
=
0.75
π
−
0.55
π
=
0.2
π
\Delta \omega = {\omega _{st}} - {\omega _p} = 0.75\pi - 0.55\pi = 0.2\pi
Δω=ωst−ωp=0.75π−0.55π=0.2π。通过查表,海明窗过渡带为
Δ
ω
=
6.6
π
/
N
\Delta \omega = 6.6\pi /N
Δω=6.6π/N
所以可知
N
=
6.6
π
Δ
ω
=
6.6
π
0.2
π
=
33
N = \frac{{6.6\pi }}{{\Delta \omega }} = \frac{{6.6\pi }}{{0.2\pi }} = 33
N=Δω6.6π=0.2π6.6π=33
由此群延迟
τ
\tau
τ:
τ
=
N
−
1
2
=
16
\tau = \frac{{N - 1}}{2} = 16
τ=2N−1=16
因此,海明窗
ω
(
n
)
\omega (n)
ω(n)为:
ω
(
n
)
=
[
0.54
−
0.46
cos
(
2
π
n
N
−
1
)
]
R
N
(
n
)
=
[
0.54
−
0.46
cos
π
n
16
]
R
33
(
n
)
\omega (n) = \left[ {0.54 - 0.46\cos (\frac{{2\pi n}}{{N - 1}})} \right]{R_N}(n){\rm{ = }}\left[ {0.54 - 0.46\cos \frac{{\pi n}}{{16}}} \right]{R_{33}}(n)
ω(n)=[0.54−0.46cos(N−12πn)]RN(n)=[0.54−0.46cos16πn]R33(n)
则线性相位FIR低通滤波器的h(n)为
接下来验证
H
(
e
j
ω
)
=
Σ
n
=
0
32
h
(
n
)
e
−
j
ω
n
H\left( {{e^{j\omega }}} \right) = \mathop \Sigma \limits_{n = 0}^{32} h(n){e^{ - j\omega n}}
H(ejω)=n=0Σ32h(n)e−jωn。此处采用MATLAB工具画出滤波器幅度响应如图。
由此可以看出阻带衰减还没有达到50dB的要求。所以考虑改选N=35的海明窗进行验证。
ω
(
n
)
=
[
0.54
−
0.46
cos
(
2
π
n
N
−
1
)
]
R
N
(
n
)
=
[
0.54
−
0.46
cos
π
n
17
]
R
35
(
n
)
\omega (n) = \left[ {0.54 - 0.46\cos (\frac{{2\pi n}}{{N - 1}})} \right]{R_N}(n){\rm{ = }}\left[ {0.54 - 0.46\cos \frac{{\pi n}}{{17}}} \right]{R_{35}}(n)
ω(n)=[0.54−0.46cos(N−12πn)]RN(n)=[0.54−0.46cos17πn]R35(n)
绘制h(n)和滤波器幅度响应、相位响应如图。
根据程序中对阻带最小衰减的判断确定,As最小为53.22dB,满足了衰减超过50dB的要求。能够有效的抑制
x
3
(
t
)
{x_3}(t)
x3(t)的频谱分量。
clc;
clear all;
fp=1100;fst=1500;%设置通带和阻带频率
Ts=0.00025; %采样间隔为0.25ms
Fs=1/Ts; %采样频率
wp=2*pi*fp/Fs;ws=2*pi*fst/Fs;%数字角频率
deltaw=ws-wp; %求过度带宽
N0=ceil(6.6*pi/deltaw);%计算海明窗长
N=N0+mod(N0+1,2);%取N为奇数
n=[0:N-1];
wd=(hamming(N))';%求海明函数
wc=(ws+wp)/2;%理想低通滤波器截至频率
hd=ideallp(wc,N);%理想低通的单位冲击响应,调用ideallp将其放于相同路径下
h=hd.*wd;
[db,mag,pha,grd,w] = freqz_m( h,[1]);%检查设计结果的各项指标(“dB”是负值),调用freqz_m将其放于相同路径下
dw=2*pi/1000;
Rp=-min(db(1:wp/dw +1));%检查通带最大衰减是否满足要求
As=-max(db(ws/dw+1:501));%检查阻带最小衰减是否满足要求
figure(1); %绘制海明窗
plot(n,wd,'linewidth',2);grid;
title('海明窗'); xlabel( 'n'); ylabel( 'w(n) ');
axis([0,N,0,1]);
figure(2);%绘制幅度响应曲线
plot(w/pi,db,'linewidth',2)
title('幅度响应(db) ' );xlabel('\omega/\pi');ylabel('20log|H(e^j^\omega)|(dB)');
axis([0,1,-80,5]);grid;
set(gca,'xtickmode','manual','xtick',[ 0,0.55,0.75,1.0]);%横坐标设置,观察感兴趣的位置,此处选取通带和阻带位置
figure(3);%绘制相位响应曲线
plot(w/pi,pha,'linewidth',2);axis([0,1,-4,4]);grid
title('相位响应');xlabel ( '\omega/\pi'); ylabel ( 'arg20log[H(e^j^\omega)]');
调用的ideallp和freqz_m函数如下(将他们放于同一工作目录下即可调用):
function hd = ideallp(wc,N);
% 理想低通滤波器的脉冲响应子程序
% hd = ideallp(wc,N)
% hd = 点0 到 N-1之间的理想脉冲响应
% wc = 截止频率(弧度)
% N = 理想滤波器的长度
%
tao = (N-1)/2; % 理想脉冲响应的对称中心位置
n = [0:(N-1)]; % 设定脉冲响应长度
m = n-tao+eps; % 加一个小数以避免零作除数
hd =sin(wc*m) ./ (pi*m); % 理想脉冲响应
function [db,mag,pha,grd,w]=freqz_m(b,a)
%求取系统的绝对幅度响应、相对的db值幅度响应、相位响应和群延时响应的函数
%db为相对振幅(dB)
%mag为绝对振幅
%pha为相位响应
%grd为群延时
%w为频率样本点向量
%b为Ha(z)分子多项式系数(对FIR而言,b=h)
%a为Hz(z)分母多项式系数(对FIR而言,a=1)
% [H,w]=freqz(b,a,1000,'whole');
%freqz显示数字滤波器频域中的图形
%[H,W] = FREQZ(B,A,N,'whole') uses N points around the whole unit circle.
[H,w] = freqz(b,a,1000,'whole');
H = (H(1:1:501))'; w = (w(1:1:501))';
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
pha = angle(H);
% pha = unwrap(angle(H));
grd = grpdelay(b,a,w);
参考文章
数字信号处理教程第四版_程佩青P210(例3.20) P545 P560(例8.1 例8.7)
链接: 数字信号处理及MATLAB实现
链接: Matlab快速入门(五)傅里叶变换