数字信号与图像处理实验二：快速傅里叶变换
实验二：快速傅里叶变换
掌握信号处理的基本思想，理解采样信号的频谱特性，加强信号采样与重建的有关基本概念的理解，深入理解线性时不变系统输出与输入的关系，了解数字信号采样率转换前后信号频谱的特征。
| 班级 | 姓名 | 实验名 | 完成时间 |
|---|---|---|---|
| 计算2014 | lx | 快速傅里叶变换 | 2022.11.14 |
1. 给定序列，并对其进行相应的分析
1.1 高斯序列
x_a(n)=\begin{cases} e^{-\frac{(n-p)^2}{q}}~~ 0\leq n\leq 15 \\ ~~ ~~ ~~0 ~~ ~~ ~~ ~~~~ ~~其他 \\ \end{cases}
观察高斯序列的时域和幅频特性，固定信号中参数p=8，改变的值，使q分别等于2， 4， 8，观察它们的时域和幅频特性，了解当q取不同值时对信号序列的时域和幅频特性的影响；固定q=8，改变p，使p分别等于8， 13， 14，观察参数p变化对信号序列的时域及幅频特性的影响，记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线，并说明频谱出现了哪些现象？
clc;clear;
n = 0:15; p1 = 8; p2 = 13; p3 = 14; q1 = 2; q2 = 4; q3 = 8;
% p = 8, q = 2, 4, 8
x1 = exp(-(n-p1).*(n-p1)/q1); x2 = exp(-(n-p1).*(n-p1)/q2); x3 = exp(-(n-p1).*(n-p1)/q3);
x1w = fft(x1); x2w = fft(x2); x3w = fft(x3);
figure(1);
subplot(3,2,1); stem(n, x1, 'fill'); ylabel('p=8, q=2'); title('时域');
subplot(3,2,2); stem(n, abs(x1w), 'fill'); title('频域');
subplot(3,2,3); stem(n, x2, 'fill'); ylabel('p=8, q=4')
subplot(3,2,4); stem(n, abs(x2w), 'fill');
subplot(3,2,5); stem(n, x3, 'fill'); ylabel('p=8, q=6')
subplot(3,2,6); stem(n, abs(x3w), 'fill');
sgtitle('高斯序列时域和幅频特性');
% q = 8, p = 8, 13, 14
x1 = exp(-(n-p1).*(n-p1)/q3); x2 = exp(-(n-p2).*(n-p2)/q3); x3 = exp(-(n-p3).*(n-p3)/q3);
x1w = fft(x1); x2w = fft(x2); x3w = fft(x3);
figure(2);
subplot(3,2,1); stem(n, x1, 'fill'); ylabel('p=8, q=8'); title('时域');
subplot(3,2,2); stem(n, abs(x1w), 'fill'); title('频域');
subplot(3,2,3); stem(n, x2, 'fill'); ylabel('p=13, q=8')
subplot(3,2,4); stem(n, abs(x2w), 'fill');
subplot(3,2,5); stem(n, x3, 'fill'); ylabel('p=14, q=8')
subplot(3,2,6); stem(n, abs(x3w), 'fill');
sgtitle('高斯序列时域和幅频特性');
由上图比较可知当保持p不变时，随q的增大时域波形展宽且变得平滑，而频域波形变陡，且频谱分量变少，不容易发生混叠，可见高斯序列中q表示时域波形的陡峭程度。
由上图比较可知当保持q不变时，随p的增大时域波形整体向右移动，而且当时域的窗口固定时在p=13和p=14有几点被泄漏了，可见p值表示时域波形峰值的位置，p值越大泄漏的越多，而频域波形随p的增大频率分量会增多，容易产生混叠。
1.2 衰减正弦序列
x_b(n)=\begin{cases} e^{-an}sin(2\pi fn)~~ 0\leq n\leq 15 \\ ~~ ~~ ~~ ~~0 ~~ ~~ ~~ ~~ ~~其他 \\ \end{cases}
观察衰减正弦序列的时城和幅频特性，a=0.1，f= 0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线。改变f，使f分别等于0.4375和0.5625，观察两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象，哪种满足采样定理。
n=0:1:15; a=0.1; f1=0.0625; f2=0.04375; f3=0.05625;
xb1=exp(-a*n).*sin(2*pi*f1*n);
figure
subplot(3, 2, 1); stem(n, xb1,'.');
title("f=0. 0625的时域特性"); xlabel('n'); ylabel('xb1 (n)'); grid on;
[H, w] = freqz(xb1, 1, [], 'whole', 1);
Hamplitude = abs (H);
subplot(3, 2, 2); plot (w, Hamplitude);
title("f=0. 0625的幅频响应"); xlabel('w/ (2*pi)'); ylabel('|H(exp(jw))|');
grid on
xb2=exp(-a*n) .* sin(2*pi*f2*n) ;
subplot(3, 2, 3); stem(n, xb2, '.');
title("f=0. 04375的时域特性"); xlabel('n'); ylabel('xb2 (n)'); grid on
[H, w] = freqz(xb2, 1, [], 'whole', 1);
Hamplitude = abs (H);
subplot(3, 2, 4);
plot(w, Hamplitude);
title("f=0. 04375的幅频响应"); xlabel('w/(2*pi)' ); ylabel('|H(exp(jw))|'); grid on
xb3=exp(-a*n) .* sin (2*pi*f3*n);
subplot (3, 2, 5);
stem(n, xb3,'.' );
title("f=0. 05625的时域特性"); xlabel('n'); ylabel('xb3(n)');
grid on
[H, w] = freqz(xb3, 1, [], 'whole', 1);
Hamplitude = abs (H);
subplot(3, 2, 6); plot (w,Hamplitude);
title("f=0. 05625的幅频响应"); xlabel('w/(2*pi)'); ylabel('|H(exp(jw))|'); grid on
当f=0.0625时，谱峰位置出现正确，但是采样频率 fs = 1/T 小于两倍的信号频率，存在在混叠现象，时域采样为一周期，不满足采样定理。
满足Nyquist定理时，f≤0.5，f=0.5625 时不满足Nyquist定理。随着f的增大，频谱的谱峰逐渐向右平移，两谱峰逐渐向中间靠拢。因为0.4375=0.5-0.0625,0.5625=0.5+0.0625，f=0.4375 和f=0.5625频谱图关于x轴对称，造成观察到的频谱完全相同，但实际上表示的意义却不相同。由于存在泄漏现象，出现了高频分量，虽然在f=0.4375时满足Nyquist定理但实际上已发生了频谱混叠。
1.3 三角波序列
x_c(n)=\begin{cases} n~~~~~~~~~~~0\leq n\leq 3 \\ 8-n~~~~4\leq n\leq 7 \\ 0~~~~~~~~~~~其他 \\ \end{cases}
1.4 2反三角波序列
x_d(n)=\begin{cases} 4-n~~~~0\leq n \leq 3 \\ n-4~~~~4\leq n \leq 7 \\ 0~~~~~~~~~~~其他 \end{cases}
观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT 分析信号序列三角波和反三角波序列的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线，说明三角波和反三角波的圆周位移关系。
close all; clear; clc;
n1 = 0:1:3; xc1 = n1+1; n2 = 4:7; xc2 = 8-n2;
xc = [xc1,xc2]; n = [n1,n2];
subplot(321)
stem(n,xc); xlabel ('n'); ylabel('xc'); title('三角序列');
n1 = 0:1:3; xd1 = 4-n1; n2 = 4:7; xd2 = n2-3;
xd = [xd1, xd2]; n = [n1,n2];
subplot(322)
stem(n,xd); xlabel('n') ; ylabel ('xd') ; title('反三角序列');
N=16;
[H1,w1] = freqz (xc,1, 256, 'whole', 1);
Hamplitude1 = abs (H1);
subplot(323)
plot (2*w1, Hamplitude1);
title('xc幅频响应'); xlabel ('w/pi'); ylabel ('IH(exp(jw)) 1'); grid on
[H2,w2] = freqz (xd,1, 256, 'whole', 1) ;
Hamplitude2 = abs (H2);
subplot(324)
plot(2*w2, Hamplitude2);
title('xd幅频响应'); xlabel ('w/pi'); ylabel('|H(exp(jw))|'); grid on
[H3, w3] = freqz(xc, 1, N, 'whole', 1);
Hamplitude3 = abs (H3) ;
subplot(325);
h3 = stem(2*w3, Hamplitude3, '*');
title ( 'xc幅频响应进行N点FFT'); xlabel ('n'); ylabel(' |H(exp(jw))|'); grid on
[H4, w4] = freqz(xd, 1, N, 'whole', 1) ;
Hamplitude4 = abs (H4);
subplot(326);
h4 = stem(2*w4, Hamplitude4, '*');
title ( 'xd幅频响应进行N点FFT'); xlabel ('n'); ylabel('|H(exp(jw))|'); grid on
由图可知，N=8时正、反三角波的序列互补，频域图形相同。因为作DFT时要先周期延拓作完后取主值部分，而正反三角波周期延拓后是相同的，只差一个相位，因此得到的频域图形也是相同的。分析三角和反三角的序列可以发现，满足圆周移位关系，三角序列经过时移4个单位后，与反三角序列相同。
2. 分析有限长正弦序列
分别用离散时间傅里叶变换(DTFT)和离散傅里叶变换(DFT)分析一个有限长正弦序列的频谱，说明两个频谱的区别。信号频率分别为10 Hz和11 Hz，采样频率为64 Hz，信号长度为32点。
clear all; clc; close all;
f1=10; f2=11; fs=64; N=32;
dt=1/fs;
T=(0:N-1)*dt;
x1=sin(2*pi*f1*T);
x2=sin(2*pi*f2*T);
x1_dft=fft(x1); x2_dft=fft(x2);
subplot(221); plot(T,abs(x1_dft/N), '-*'); xlabel('频率/HZ'); ylabel('振幅'); title('f=10HZ DFT变换');
subplot(222); plot(T,abs(x2_dft/N), '-*'); xlabel('频率/HZ'); ylabel('振幅'); title('f=11HZ DFT变换');
[h1,w1]=freqz(x1,1); [h2,w2]=freqz(x2,7);
subplot(223); plot(w1,abs(h1)); xlabel('频率/HZ'); ylabel('振幅'); title('f=10HZ DTFT变换');
subplot(224); plot(w2,abs(h2)); xlabel('频率/HZ'); ylabel('振幅'); title('f=11HZ DTFT变换');
观察两个序列的频谱图，发现两张频谱图都是对称的，但是f=11HZ的频谱图明显有很多低频分量，出现了信号泄露。两个序列的离散时间傅里叶变换相同。
为什么DFT是DTFT的等间隔采样？由于DFT是DTFT的取样值，其相邻两个频率样本点的间距2πN，所以如果我们增加数据的长度N，使得到的DFT谱线就更加精细，其包络就越接近DTFT的果，这样就可以利用DFT计算DTFT。如果没有更多的数据，可以通过补零来增加数据长度。
实验要求：
- 熟悉MATLAB的信号处理函数
- 熟练离散信号时域和频域图形的绘制
- 熟练绘制信号的幅频响应和相频响应
- 完成频谱与实例分析
附：相关MATLAB函数
y =exp( X)；对向量X的各元素作指数计算，结果为一向量。
conj(X)；对向量X的各元素作复共辄运算，即将虚部改变符号。
real(X)；对向量X的各元素取其实部。
v =randn ( size(X))；生成同X具有相同维数的正态(高斯)分布的随机矩阵，通常用于生成高斯白噪声。
fft(X，N)；计算X的N点FFT，如果X的长度小于N，则将在X后补零。反之，如果X的长度大于N，则将对X进行截取。
ifft(X)；计算X的N点 IFFT。
fftshift(Y)；如果Y为向量，该命令将Y分成左右两部分并交换位置。
xcorr(A ，B)；计算A、B两个序列的互相关。
n ( size(X))；生成同X具有相同维数的正态(高斯)分布的随机矩阵，通常用于生成高斯白噪声。fft(X，N)；计算X的N点FFT，如果X的长度小于N，则将在X后补零。反之，如果X的长度大于N，则将对X进行截取。
ifft(X)；计算X的N点 IFFT。
fftshift(Y)；如果Y为向量，该命令将Y分成左右两部分并交换位置。
xcorr(A ，B)；计算A、B两个序列的互相关。