实验一 模拟信号频谱分析
1.实验目的
- 学会应用DFT对模拟信号进行频谱分析的方法;
- 通过应用DFT分析各种模拟信号的频谱,加深对DFT的理解;
- 熟悉MATLAB的基本操作,以及一些基本函数的使用,为以后的实验奠定基础。
3.实验内容
⑵ 假设一实际测得的一段信号的长度为0.4秒,其表达式为:
其中。试确定一合适抽样频率,利用MATLAB的fft函数分析计算信号的频谱。
解:信号的最高频率,抽样频率,取抽样频率;最低的频率分辨率为,最少的信号样点数为。
程序:????????????????????????????????????????????结果:
N=30; %数据的长度
L=512; %DFT的点数
f1=100;
f2=110;
fs=300; %抽样频率
T=1/fs; %抽样间隔
t=(0:N-1)*T;
x=cos(2*pi*f1*t)+0.75*cos(2*pi*f2*t);
y=fft(x,L);
mag=abs(y);
f=(0:length(y)-1)‘*fs/length(y);
plot(f(1:L/2),mag(1:L/2));
xlabel(‘频率(Hz)‘)
ylabel(‘幅度谱‘)
?
⑶ 观察并分析采用不同抽样频率时,对信号的频谱影响。
因为:
?
其幅度频谱为
- 以,对其进行采样得到;
程序: N=256; %数据的长度
fs=5000; %抽样频率
Ts=1/fs; %抽样间隔
t1=(0:N-1)*Ts;
t2=(-N+1:0)*Ts;
x=exp(1000*t2)+exp(-1000*t1);
y=fft(x);
mag=Ts*abs(y);
Ws=2*pi*fs;
w=(0:length(y)-1)‘*Ws/length(y);
w1=w(1:length(y)/2);
plot(w1,mag(1:length(y)/2),‘r‘);
xlabel(‘频率(弧度/秒)‘);
ylabel(‘幅度谱‘);
z=[‘N=‘ num2str(N) ‘ fs=‘ num2str(fs) ‘的结果‘];
legend(z);
title(‘exp(-t)的幅度谱‘)
结果:
- 以,对其进行采样得到。
程序:N=256; %数据的长度
fs=1000; %抽样频率
Ts=1/fs; %抽样间隔
t1=(0:N-1)*Ts;
t2=(-N+1:0)*Ts;
x=exp(1000*t2)+exp(-1000*t1);
y=fft(x);
mag=Ts*abs(y);
Ws=2*pi*fs;
w=(0:length(y)-1)‘*Ws/length(y);
w1=w(1:length(y)/2);
plot(w1,mag(1:length(y)/2),‘r‘);
xlabel(‘频率(弧度/秒)‘);
ylabel(‘幅度谱‘);
z=[‘N=‘ num2str(N) ‘ fs=‘ num2str(fs) ‘的结果‘];
title(‘exp(-t)的幅度谱‘)
?
结果:
?
实验二 离散信号频谱分析
1.实验目的
- plot(w1,X(1:length(y)/2),‘-‘,w1,mag(1:length(y)/2),‘r‘);理解DFS、IDFS的原理和基本性质;
- 掌握应用FFT对离散信号进行频谱分析的方法;
- 通过应用FFT分析各种离散信号的频谱,学会在实际中正确应用FFT。
3.实验内容
⑴ 理解运行以上例题程序,改变有关参数,进一步观察结果的变化,并加以分析说明。
⑵ 已知序列: ,试确定一合适样本数,利用MATLAB的fft函数分析计算信号的频谱。
10点DFT程序:
n = [0:1:9];
x = cos(0.82*pi*n)+2*sin(0.43*pi*n);
subplot(2,1,1);
stem(n,x,‘.‘);
title(‘x(n),0 <= n <= 9‘);
xlabel(‘n‘);
ylabel(‘x(n)‘);
axis([0,10,-2.5,2.5]);
Xk = fft(x);
magXk = abs(Xk(1:1:6));
k1 = 0:1:5;
w1 = 2*pi/10*k1;
subplot(2,1,2);
stem(w1/pi,magXk,‘.‘);
title(‘DTFT采样‘);
xlabel(‘频谱‘);
ylabel(‘|X(k)|‘);
axis([0,1,0,10]);
?
结果如图:
?
?
?
?
?
⑶ 利用FFT计算信号的频谱。
a=0.5;
Nsum=40;
n=0:Nsum-1;
for i=1:4
switch i
case 1
N=5;
case 2;
N=10;
case 3
N=20;
case 4
N=40;
end
k=0:N-1;
zk=exp(j*2*pi*k/N);
Xk=zk./(zk-a);
xnN=real(ifft(Xk,N));
xn=xnN‘*ones(1,Nsum/N);
xn=(xn(:))‘;
gsptext=strcat(‘22‘,int2str(i));
subplot(gsptext)
stem(n,xn,‘.‘);axis([0,40,-0.1,1.5]);
xtext=strcat(‘IDFT N=‘,int2str(N));
xlabel(‘n‘);ylabel(‘x(n)‘);title(xtext);
hold on
plot(n,zeros(1,length(n)));
hold off
End
结果
如图:
?
实验三 IIR数字滤波器的设计
1.实验目的
- 熟悉IIR数字滤波器的设计方法;
- 掌握IIR数字滤波器的MATLAB实现程序;
- 进一步加深对数字滤波器的常用指标和设计过程的理解;
- 通过观察对实际心电图信号的滤波作用,获得数字滤波工程应用的感性知识。
3.实验内容
⑴ 运行以上例题程序,加深对上述几种IIR数字滤波器的常用指标和设计过程及MATLAB实现程序的理解。
⑵ 人体心电图信号在测量过程中往往受到工业高频干扰,所以必须经过低通滤波处理后,才能得到判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n),其中存在高频干扰。试设计一个低通滤波器,滤除其中的干扰成分。设计指标:通带边界频率ωp=40Hz,通带波纹小于0.5dB,阻带边界频率ωs=50Hz,阻带衰减大于40dB。采样频率fs=200Hz。
- 设计一个切比雪夫Ⅰ型低通滤波器并绘出它的幅频响应曲线;
Fs=200; %采样频率
wp=40*2/Fs; %通带边界频率
ws=50*2/Fs; %阻带边界频率
Rp=0.5; %通带波纹小于0.5dB
Rs=40;???????????? %阻带衰减大于40dB
Nn=128;
[N,Wn]=cheb1ord(wp,ws,Rp,Rs)
[b,a]=cheby1(N,Rp,Wn);
filt(b/a(1),a/a(1),1/Fs)
figure(1)
freqz(b,a,Nn,Fs)
figure(2)
[H,F]=freqz(b,a,Nn,Fs);
magH=abs(H);
plot(F,magH);
xlabel(‘Frequency (Hz)‘);
ylabel(‘Magnidute‘);
title(‘Chebyshev I type digital bandpass filter‘)
Grid
结果:N = 8
Wn = 0.4000
Transfer function:
0.0003471 + 0.002777 z^-1 + 0.009719 z^-2 + 0.01944 z^-3 + 0.0243 z^-4 + 0.01944 z^-5
+ 0.009719 z^-6 + 0.002777 z^-7 + 0.0003471 z^-8
--------------------------------------------------------------------------------------
1 - 3.866 z^-1 + 8.262 z^-2 - 11.69 z^-3 + 11.78 z^-4 - 8.544 z^-5 + 4.356 z^-6
- 1.434 z^-7 + 0.2381 z^-8
Sampling time: 0.005
响应曲线如图:
?
?
- 用设计的滤波器对原数据序列进行滤波;
- 绘制滤波以后的信号波形,并与原信号进行比较。
x(n)=[ -4,-2, 0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2, 6, 12, 8, 0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4, 8, 12, 12, 10, 6, 6, 4, 0, 0, 0, 0, 0,-2,-2, 0, 0,-2,-2,-2,-2, 0 ]
?
Fs=200; %采样频率
wp=40*2/Fs; %通带边界频率
ws=50*2/Fs; %阻带边界频率
Rp=0.5; %通带波纹小于0.5dB
Rs=40;???????????? %阻带衰减大于40dB
Nn=128;
n =[0:1:49];
x=[ -4,-2, 0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2, 6, 12, 8, 0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4, 8, 12, 12, 10, 6, 6, 4, 0, 0, 0, 0, 0,-2,-2, 0, 0,-2,-2,-2,-2, 0 ]
[N,Wn]=cheb1ord(wp,ws,Rp,Rs)
[b,a]=cheby1(N,Rp,Wn);
filt(b/a(1),a/a(1),1/Fs)
figure(1)
y=filter(b,a,x)
subplot(221);
plot(n,y)
axis([0 50 -100 20])
title(‘滤波后波形‘)
grid
subplot(222);
plot(n,x)
axis([0 50 -100 20])
title(‘原波形‘)
grid
?
结果如图: