华工信号与系统matlab实验[实验三--利用DFT分析模拟信号频谱]

信号与系统[实验三–利用DFT分析模拟信号频谱]
实验题目

1. 利用FFT分析信号 的频谱。
   (1) 确定DFT计算的各参数（抽样间隔，截短长度，频谱分辨率等）；
   (2) 比较理论值与计算值，分析误差原因，提出改善误差的措施。
2. 分析周期信号的频谱时，如果分析长度不为整周期，利用fft函数计算并绘出其频谱，总结对周期信号进行频谱分析时，如何选取信号的分析长度。
3. 假设一实际测得的一段信号的长度为0.4s，其表达式为
   x(t)=cos(2pf1t)+0.75 cos(2pf2t)
   其中f1=100Hz，f2=110Hz。当利用FFT近似分析该信号的频谱时，需要对信号进行时域抽样。试确定一合适抽样频率fsam，利用DFT分析信号x(t)的频谱。
   若在信号截短时使用Hamming窗，由实验确定能够分辨最小谱峰间隔Df和信号长度Tp的关系。若采用不同参数的Kaiser窗，重新确定能够分辨最小谱峰间隔Df和信号长度Tp 的关系。
4. 产生一个淹没在噪声中的信号x(t)，例如由50Hz和120Hz的正弦信号以及一个零均值的随机噪声叠加而成。确定抽样间隔和信号截短长度，分析信号的频谱，指出50Hz和120Hz的正弦成分对应的谱峰位置，详细写出检测信号的步骤和原理。

第一题代码:

fsam=50;
Tp=6;
N=512; 
T=1/fsam;
t=0:T:Tp; 
x=exp(-2*t); 
X=T*fft(x,N);        %频域幅度  去归一化,乘以采样间隔T

subplot(3,1,1);
plot(t,x);
xlabel('t');
ylabel('|X(jw)|')
title('时域波形');

w=(-N/2:N/2-1)*(2*pi/N)*fsam;    %频域频率  去归一化,除以T
y=1./(j*w+2);

subplot(3,1,2);
plot(w,abs(fftshift(X)),w,abs(y),'r-.');
title('幅度谱');
xlabel('w');
ylabel('y=|1/(j*w+1)|')
legend('理论值','计算值');
axis([-10,10,0,1.4])

subplot(3,1,3);
error = abs(abs(fftshift(X))-abs(y));
plot(w,error)
title('理论与计算的误差');

1)抽样间隔 0.02
截短长度 6
频谱分辨率 50/512
2)由图3-1-1可以看出理论值与实际值相差不多.误差原因是采样间隔不够小.
改进:减小采样间隔

第二题代码:

close all;
%题二 ---连续周期函数
T0=1;       %w0 = 2*pi,  T0 = 1  [w0 = 最大公约数]
T1 = 1.5;   %分析长度不为整周期

w0 = 2*pi ; %最大公约数
p = 9; %最高谐次  110/10
N = 99;  %N>2*p+1
T = T0/N;     %周期T0、FFT的点数N、抽样间隔T
T2 = T1/N;

t = 0:T:T0;
t1 = 0:T2:T1;

x=cos(2*pi*5*t)+2*sin(2*pi*9*t);  %周期信号 
x1=cos(2*pi*5*t1)+2*sin(2*pi*9*t1); 


figure
subplot(2,1,1)
Xm=fft(x,N) /N;  %频域幅度  得到的是傅里叶级数
w=(-(N-1)/2:(N-1)/2)*(2*pi/T0);   %T0为1个周期,即显示-pi  ~  pi
stem(w,abs(fftshift(Xm)));
xlabel('w (rad)');
ylabel('magnitude');  
title('分析长度为整周期'); 


subplot(2,1,2)
Xm1=fft(x1,N)/N;  %利用FFT计算其频谱
w=(-(N-1)/2:(N-1)/2)*(2*pi/T1);
stem(w,abs(fftshift(Xm1)));
xlabel('w (rad)');
ylabel('magnitude');  
title('分析长度为非整周期'); 

第三题代码:

close all;
clear all;

%连续非周期
fsam = 220; %抽样频率,cos经过傅里叶变换相当于频率搬移,f2搬移到-110  +110,为脉冲信号
            %fsam = 2*fm
T0 = 0.4;
f01 = 100;
f02 = 110;

w0 = 10*pi ; %最大公约数
p = 11; %最高谐次  110/10

%T=1/fsam;     %周期T0、FFT的点数N、抽样间隔T  
T = 0.004;
t=0:T:(T0-T);
N = length(t);
x=cos(2*pi*f01*t)+ 0.75*cos(2*pi*f02*t);

%问1-------------------------------------
figure
subplot(2,1,1)
plot(t,x)
xlabel('t');
title('x(t)');

subplot(2,1,2)
Xm=fft(x,N)/N;  %傅里叶级数
w=(-N/2:N/2-1)*(2*pi/N)/T0;
f = ((-N/2 + 1):N/2)/T0;
stem(f,abs(fftshift(Xm)));
xlabel('f(hz)');
ylabel('magnitude');  

%问2-------------------------------------
%在信号截短时使用Hamming窗，由实验确定能够分辨最小谱峰间隔Df和信号长度Tp的关系

%实验验证DF与Tp无关,无论Tp取何值都是相差10HZ,前提是Tp至少有一个周期长度
T2 = 0.004;
T02 = 0.4;
t2 = 0:T2:(T02-T2);
N2 = length(t2)
k2 = 1:N2;

Ham2 = 1/2.*(1-cos(2*pi* k2 /(N2-1)));
x2 = (cos(2*pi*f01*t2)+ 0.75*cos(2*pi*f02*t2)).*Ham2;

figure
subplot(2,1,1)
plot(t2,x2,t2,Ham2,'r-')
legend('x(t)','Ham')
xlabel('t');
title('x(t) and Hamming窗');

subplot(2,1,2)
Xm2=fft(x2,N2)/N2;  %傅里叶级数

f2 = ((-N2/2 + 1):N2/2)/T02;
stem(f2,abs(fftshift(Xm2)));
xlabel('f(hz)');
ylabel('magnitude');

%问3-------------------------------------
%在信号截短时使用Kaiser窗，由实验确定能够分辨最小谱峰间隔Df和信号长度Tp的关系
%实验验证DF与Tp无关,无论Tp取何值都是相差10HZ,前提是Tp至少有一个周期长度
T3 = 0.004;
T03 = 0.4;
t3 = 0:T3:(T03-T3);
N3 = length(t3);
k3 = 1:N3;

beta = 10;
Kai=kaiser(N3,beta);

x3 = (cos((2*pi*f01).*t3)+ 0.75*cos((2*pi*f02).*t3)).*Kai';

figure
subplot(2,1,1)
plot(t3,x3,t3,Kai)

legend('x(t)','Kaiser')
xlabel('t');
title('x(t) and Kaiser窗');
ylabel('x(t)')

subplot(2,1,2)
Xm3=fft(x3,N3)/N3;  %傅里叶级数

f3 = ((-N3/2 + 1):N3/2)/T03;
stem(f3,abs(fftshift(Xm3)));
xlabel('f(hz)');
ylabel('magnitude');

理论上fsam应取220,可是选择220分析的时候,采样间隔为无限循环小数,导致fft后会丢失110Hz的脉冲,所以选择为fsam = 250
无论是Hamming窗还是Kaiser窗,通过实验改变分析长度Tp,得到的结果都是,两个谱峰间的间隔都是10Hz.故能够分辨最小谱峰间隔Df与信号长度Tp 无关.

第四题代码:

% 题4
close all;
clear all;

fsam = 250;
T0 = 1;
T = 1/fsam;
t = 0:T:(T0-T);
N=length(t); 
f1 = 50;
f2 = 120;

x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + randn(size(t));

subplot(211)
plot(t,x)
xlabel('t');
ylabel('x(t)');
title('含随机干扰的信号x(t)');


X=T*fft(x,N);        %频域幅度  去归一化,乘以采样间隔T
f=(-(N-1)/2:(N-1)/2)/T0;    %频域频率  去归一化,除以T

subplot(212)
plot(f,abs(fftshift(X)));
title('幅度谱');
xlabel('f(Hz)');
ylabel('magnitude');

由于f1 = 50Hz, f2 = 120Hz,所以fsam = 240Hz
可是如果取fsam = 240Hz,采样间隔为无限循环小数,所以选取为250Hz,
则
抽样间隔 = 1/250 =0.004 s
信号截短长度 1s
50Hz和120Hz的频谱可在图3-4-1看出
检测信号的步骤和原理:
由于x(t)是混频相加的函数,可通过傅里叶变换得到各频率上的脉冲信号,即可检测出.
本题多了干扰函数,但是干扰函数的幅度较小,不会对在50Hz和120Hz的脉冲信号有太大的影响,所以可以分辨出.