实验四：用 FFT 进行谱分析

实验内容

绘制出离散幅度谱

%% 1. 对信号 x1(n) 做 8 点和 16 点的 FFT
x = [1, 1, 1, 1, 0, 0, 0, 0];
N = 8;
Plot_waveforms_spectra(x, N, 'x1 的波形', 'x1 的 8 点FFT');
N = 16;
Plot_waveforms_spectra(x, N, 'x1 的波形', 'x1 的 16 点FFT');
%% 2. 对信号 x2(n)、x3(n) 做 8 点和 16 点的 FFT
x = [1, 2, 3, 4, 4, 3, 2, 1];
N = 8;
Plot_waveforms_spectra(x, N, 'x2 的波形', 'x2 的 8 点FFT');
N = 16;
Plot_waveforms_spectra(x, N, 'x2 的波形', 'x2 的 16 点FFT');
x = [4, 3, 2, 1, 1, 2, 3, 4];
N = 8;
Plot_waveforms_spectra(x, N, 'x3 的波形', 'x3 的 8 点FFT');
N = 16;
Plot_waveforms_spectra(x, N, 'x3 的波形', 'x3 的 16 点FFT');
%% 3. 对信号 x4(n)、x5(n) 做 8 点和 16 点的 FFT
n = 0:7;
x = cos(pi / 4 .* n);
N = 8;
Plot_waveforms_spectra(x, N, 'x4 的波形', 'x4 的 8 点FFT');
n = 0:15;
x = cos(pi / 4 .* n);
N = 16;
Plot_waveforms_spectra(x, N, 'x4 的波形', 'x4 的 16 点FFT');
n = 0:7;
x = sin(pi / 8 .* n);
N = 8;
Plot_waveforms_spectra(x, N, 'x5 的波形', 'x5 的 8 点FFT');
n = 0:15;
x = sin(pi / 8 .* n);
N = 16;
Plot_waveforms_spectra(x, N, 'x5 的波形', 'x5 的 16 点FFT');
%% 4. 对信号 x6(n) 以 fs=64（Hz）采样后做 16、32、64 点的 FFT
N = 16;
n = 0:N-1;
fs = 64;
T = 1 / fs;
% t = n * T;
% xa = cos(8 * pi .* t) + cos(16 * pi .* t) + cos(20 * pi .* t);
xn = cos(8 * pi * n * T) + cos(16 * pi * n * T) + cos(20 * pi * n * T);
Plot_waveforms_spectra(xn, N, 'x6 的波形', 'x6 的 16 点FFT');
N = 32;
n = 0:N-1;
fs = 64;
T = 1 / fs;
% t = n * T;
% xa = cos(8 * pi .* t) + cos(16 * pi .* t) + cos(20 * pi .* t);
xn = cos(8 * pi .* n * T) + cos(16 * pi .* n * T) + cos(20 * pi .* n * T);
Plot_waveforms_spectra(xn, N, 'x6 的波形', 'x6 的 32 点FFT');
N = 64;
n = 0:N-1;
fs = 64;
T = 1 / fs;
% t = n * T;
% xa = cos(8 * pi .* t) + cos(16 * pi .* t) + cos(20 * pi .* t);
xn = cos(8 * pi .* n * T) + cos(16 * pi .* n * T) + cos(20 * pi .* n * T);
Plot_waveforms_spectra(xn, N, 'x6 的波形', 'x6 的 64 点FFT');
5、编写 matlab M 文件，读取 motherland.wav 数据，分析第 8000 至 8199 共 200 个采样点的频谱（提示是傅里叶变换）。方法：对这 200 个点数据做 N=512 的 DFT（采用 FFT 实现）。要求：画出其在[0,2π)的连续幅度谱和相位谱图。
[xn, fs] = audioread('motherland.wav'); % [音频数据, 采样频率]
xn = xn(8000:8199);
N = 512;
h = figure;
plot(8001:8000+length(xn), xn);
saveas(h, '原信号采样序列', 'svg');
Plot_waveforms_spectra(xn, N, '音频数据的波形', '音频数据的 512 点FFT');
xk = fft(xn, N);
h = figure;
subplot(2, 1, 1);
plot(linspace(0, 2, N), abs(xk)); % [0, 2pi]
xlabel('\omega (\pi rad/s)');
ylabel('|X(e^{j \omega})|');
title('幅频响应曲线 |X(e^{j \omega})|');
subplot(2, 1, 2);
plot(linspace(0, 2, N), angle(xk));
xlabel('\omega (\pi rad/s)');
ylabel('\phi(\omega)');
title('相频响应曲线 \phi(\omega)');
saveas(h, '频率响应曲线', 'svg');
6、编写 Matlab 程序，分析 lena 图像的二维频谱，调用 fft2 和 fftshift 函数实现。要求：显示时域原图、二维幅度谱图。
clc, clear;
% h = figure;
A = imread('lena.bmp'); % 读原图
% subplot(131);
h1 = figure;
imshow(A); title('原图');
fftI = fft2(A); % 二维离散傅里叶变换
A1 = abs(fftI); % 取模值
saveas(h1, 'lena原图', 'svg');
h2 = figure;
% subplot(132);
% 把幅度限定在[0,255]
B1 = (A1-min(min(A1)))/(max(max(A1))-min(min(A1)))*255;
imshow(B1); title('二维幅度谱图');
B = fftshift(B1);
saveas(h2, 'lena二维幅度谱', 'svg');
h3 = figure;
% subplot(133);
imshow(B); title('移到中心位置的二维频谱图');
saveas(h3, 'lena中心位置的二维频谱', 'svg');

思考题

1. 在 N=8 和 N=16 两种情况下，的幅频特性会相同吗？为什么？
% N=8 相同，因为 x2 = {1, 2, 3, 4, 4, 3, 2, 1}, x3 = {4, 3, 2, 1, 1, 2, 3, 4}
% x2, x3 可以通过循环时移性质互相转换
% N=16 不相同，因为
% x2 = {1, 2, 3, 4, 4, 3, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0}
% x3 = {4, 3, 2, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0}
% x2, x3 不可以通过循环时移性质互相转换
2. 如果周期信号的周期预先不知道，如何用 FFT 进行分析？
% 先截取M点进行DFT,再将截取长度扩大1倍截取,比较结果,
% 如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,
% 如果不满足误差要求就继续将截取长度加倍,重复比较,直到结果满足要求
3. 序列 x=[1,1,2,2,3,3,2,2,1,1]。
（1）对 x 进行 2 选 1 的抽取，得到序列 x1=[1,2,3,2,1]；
（2）对 x 进行 0 值内插，得到序列 x2=[1,0,1,0,2,0,2,0,3,0,3,0,2,0,2,0,1,0,1]。
（3）试使用函数 fft 分别画出 x、x1 和 x2 在[0,2π)的连续幅度谱图（提示是序列傅里叶变换的幅度谱）。
（4）写出 x1 和 x2 与 x 频谱关系的数学表达式，并解释 x1 和 x2 与 x 的幅度频谱的变化。
说明：
（1）（2）（3）如下；（4）略
N = 256;
x = [1, 1, 2, 2, 3, 3, 2, 2, 1, 1];
xk = fft(x, N);
h = figure;
plot(linspace(0, 2, N), abs(xk)); % [0, 2pi]
xlabel('\pi \omega');
ylabel('|X(e^{j \omega})|');
title('连续幅度谱图-原序列');
saveas(h, '思考题3-连续幅度谱图-原序列', 'svg');
x1 = [1, 2, 3, 2, 1];
xk = fft(x1, N);
h = figure;
plot(linspace(0, 2, N), abs(xk)); % [0, 2pi]
xlabel('\pi \omega');
ylabel('|X(e^{j \omega})|');
title('连续幅度谱图-抽取');
saveas(h, '思考题3-连续幅度谱图-抽取', 'svg');
x2 = [1, 0, 1, 0, 2, 0, 2, 0, 3, 0, 3, 0, 2, 0, 2, 0, 1, 0, 1];
xk = fft(x2, N);
h = figure;
plot(linspace(0, 2, N), abs(xk)); % [0, 2pi]
xlabel('\pi \omega');
ylabel('|X(e^{j \omega})|');
title('连续幅度谱图-内插');
saveas(h, '思考题3-连续幅度谱图-内插', 'svg');
function Plot_waveforms_spectra(x, N, title1, title2)
xk = fft(x, N);
h = figure;
subplot(211); stem(0:length(x)-1, x, '.'); title(title1);
subplot(212); stem(0:N-1, abs(xk), '.'); title(title2);
saveas(h, title2, 'svg');
end