%实验4：DFT频谱分析
%设计计算机程序，产生序列并计算序列的FFT和IFFT，绘制其幅频特性和相频特性曲线；
%模拟产生离散系统的输入序列和单位脉冲响应，利用FFT和IFFT算法计算系统的输出响应，分析FFT的计算长度对系统输出响应的影响；
%模拟产生连续时间信号，选取适当的采样频率对其采样，并用FFT算法计算其频谱，分析信号的观测时间长度、FFT的计算长度对信号频谱计算结果的影响。
程序直接复制粘贴跑，程序结果输出截图：第一个要求，略...，第二个要求，L>=N-M+1,这样才能用快速卷积计算线性卷积，主要是第三个要求，分析Sa函数，增加N，提高谱分辨率，增加Tp，减少截断误差。
%1.
xn = [1,1,1,1,1,1];%模拟输入信号
N = 6;
Xk16 = fft(xn,16);
i16 = 0:15;
Xk32 = fft(xn,32);
i32 = 0:31;
figure(1)
subplot(421);stem(i16,Xk16);title("16点DFT的幅频特性图");
subplot(422);stem(i32,Xk32);title("32点DFT的幅频特性图");
subplot(423);stem(i16,angle(Xk16));title("16点DFT的相频特性图");
subplot(424);stem(i32,angle(Xk32));title("16点DFT的相频特性图");
Xn16 = ifft(Xk16);
Xn32 = ifft(Xk32);
subplot(425);stem(i16,Xn16);title("16点IDFT的幅频特性图");
subplot(426);stem(i16,Xn16);title("32点IDFT的幅频特性图");
subplot(427);stem(i16,angle(Xn16));title("16点IDFT的相频特性图");
subplot(428);stem(i32,angle(Xn32));title("32点IDFT的相频特性图");
%2.
M = 5;
n = 0:M-1;
yn = cos(pi*n/10) + cos(2*pi*n/5);%模拟系统函数
ync = conv(xn,yn);
nc = 0:length(ync)-1;
figure(2);
subplot(421);stem(nc,ync);title("卷积法求输入序列的响应");
j = 1;
for i = -3:1:3
L = M + N -i;
Xk = fft(xn,L);;
Yk = fft(yn,L);
XYk = Xk.*Yk;
XYik = ifft(XYk);
h = 0:L-1;
j = j + 1;
subplot(4,2,j);stem(h,XYik);title("");title(sprintf('%d点快速卷积',L));
end
%结论：L-M+1点卷积结果不发生重叠
%3
t=0:pi/20:16;
Sa=sin(pi.*(t-4))./(pi*(t-4));
figure(3);
subplot(511);plot(t,Sa);grid on;
Tp1 = 8;
n1 = 0:pi/20:Tp1;
f1 = sin(pi.*(n1-4))./(pi*(n1-4));
Tp2 = 10;
n2 = 0:pi/20:Tp2;
f2 = sin(pi.*(n2-4))./(pi*(n2-4));
N1 = 12000;
N2 = 1200;
%观察时间长度的影响
F1 = fft(f1,N1);
F2 = fft(f2,N1);
k1 = 0:N1-1;
subplot(512);
plot(k1,fftshift(abs(F1))./max(abs(F1)));xlabel('k');ylabel('|F(k)|');
title('不同观测时间长度下FFT的幅频特性');
legend('Tp=8s')
subplot(513);
plot(k1,fftshift(abs(F2))./max(abs(F2)));xlabel('k');ylabel('|F(k)|');
title('不同观测时间长度下FFT的幅频特性');
legend('Tp=10s')
%计算长度的影响
F1 = fft(Sa,N1);
F2 = fft(Sa,N2);
k1 = 0:N1-1;
k2 = 0:N2-1;
subplot(514);
plot(k1,fftshift(abs(F1))./max(abs(F1)));xlabel('k');ylabel('|F(k)|');
title('不同计算长度下FFT的幅频特性');
legend('N = 12000')
subplot(515);
plot(k2,fftshift(abs(F2))./max(abs(F2)));xlabel('k');ylabel('|F(k)|');
title('不同计算长度下FFT的幅频特性');
legend('N = 1200')
