实验五：时域采样与频域采样

1、时域采样理论的验证

Tp = 64/1000; % s

A = 444.128;

a = 50 * 1.414 * pi;

w = 50 * 1.414 * pi;

fs = 2048;

T = 1 / fs;

t = 0:T:0.03;

xa = A .* exp(-a * t) .* sin(w * t);

h = figure;

subplot(2, 1, 1);

plot(t, xa);

title('波形图');

xlabel('t');

ylabel('xa(t)');

subplot(2, 1, 2);

N = 1024;

temp = fs / N * 4;

xk = fft(xa, N) * T;

plot(linspace(0, 500, N/temp), abs(xk(1:N/temp)));

title('幅频特性曲线');

xlabel('f/Hz');

ylabel('|xa(jf)|');

saveas(h, '幅频特性曲线', 'svg');

%% 1000Hz

h = figure;

fs = 1000; % Hz

T = 1 / fs; % s

N = Tp * fs; % 64

% disp(N);

N = 64;

n = 0: N - 1;

xan = A .* exp(-a * n * T) .* sin(w * n * T);

xk = T * fft(xan, N);

subplot(3, 2, 1);

stem(n, xan, '.');

axis([0 80 -100 200]);

title('F_s=1000Hz');

xlabel('n');

ylabel('xa(n)');

subplot(3, 2, 2);

k = 0:N - 1;

wk = fs / N * k;

plot(wk, abs(xk));

axis([0 fs 0 1]);

title('T*FT[xa(nT)], Fs=1000Hz');

xlabel('f(Hz)');

ylabel('幅度');

%% 200 Hz

fs = 200;

T = 1 / fs;

N = Tp * fs; % 12.8 -> 16

% disp(N);

N = 16;

n = 0: N - 1;

t = 0:0.0001:0.05;

xa = A .* exp(-a * t) .* sin(w * t);

xan = A .* exp(-a * n * T) .* sin(w * n * T);

xk = T * fft(xan, N);

subplot(3, 2, 3);

stem(n, xan, '.');

axis([0 80 -100 200]);

title('F_s=200Hz');

xlabel('n');

ylabel('xa(n)');

subplot(3, 2, 4);

k = 0:N - 1;

wk = fs / N * k;

plot(wk, abs(xk));

axis([0 fs 0 1]);

title('T*FT[xa(nT)], Fs=200Hz');

xlabel('f(Hz)');

ylabel('幅度');

%% 300 Hz

fs = 300;

T = 1 / fs;

N = Tp * fs; % 19.2 -> 32

% disp(N);

N = 32;

n = 0: N - 1;

t = 0:0.0001:0.05;

xa = A .* exp(-a * t) .* sin(w * t);

xan = A .* exp(-a * n * T) .* sin(w * n * T);

xk = T * fft(xan, N);

subplot(3, 2, 5);

stem(n, xan, '.');

axis([0 80 -100 200]);

title('F_s=300Hz');

xlabel('n');

ylabel('xa(n)');

subplot(3, 2, 6);

k = 0:N - 1;

wk = fs / N * k;

plot(wk, abs(xk));

axis([0 fs 0 1]);

title('T*FT[xa(nT)], Fs=300Hz');

xlabel('f(Hz)');

ylabel('幅度');

saveas(h, '时域采样', 'svg');

2、频域采样理论的验证

① 直接调用 MATLAB 函数 fft 计算就得到在在的 32 点频率域采样

② 抽 取的 偶 数 点 即 可 得 到在的 16 点 频 率 域 采 样， 即。

③当然也可以按照频域采样理论，先将信号 x(n) 以 16 为周期进行周期延拓，取其主值区（16 点），再对其进行 16 点 DFT(FFT),得到的就是在的 16 点频率域采样。

h = figure;

xn = [1:14 13:-1:1 0];

xk = fft(xn, 1024);

xk32 = fft(xn, 32);

xn32 = ifft(xk32, 32);

xk16 = xk32(1:2:end);

xn16 = ifft(xk16, 16);

% 1

subplot(3, 2, 1);

plot(linspace(0, 1, 512), abs(xk(1:512)));

xlabel('ω/π');

ylabel('|X(exp(jω)|');

title('FT[x(n)]');

subplot(3, 2, 2);

stem(xn, '.');

axis([1 31 0 20]);

xlabel('n');

ylabel('x(n)');

title('三角波序列[x(n)]');

% 2

subplot(3, 2, 3);

k = 0:15;

stem(k, abs(xk16), '.');

axis([0 8 0 200]);

xlabel('k');

ylabel('|X16(k)|')

title('16点频域采样');

subplot(3, 2, 4);

stem(abs(xn16), '.');

axis([1 31 0 20]);

xlabel('n');

ylabel('x16(n)');

title('16点IDFT[X_{16}(k)]');

% 3

subplot(3, 2, 5);

k = 0:31;

stem(k, abs(xk32), '.');

axis([0 16 0 200]);

xlabel('k');

ylabel('|X32(k)|')

title('32点频域采样');

subplot(3, 2, 6);

stem(abs(xn32), '.');

axis([1 31 0 20]);

xlabel('n');

ylabel('x32(n)');

title('32点IDFT[X_{32}(k)]');

saveas(h, '频域采样', 'svg');

3、语音信号整数倍抽取后频谱的分析

[x, fs] = audioread('motherland.wav'); % [音频数据, 采样频率]

%% （1）

xn = x(1000:1:2999);

N = 2048;

xk = fft(xn, N);

h = figure;

plot(8001:8000+length(xn), xn);

title('原信号波形');

saveas(h, '原信号波形', 'svg');

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');

%% （2）

xn = x(1000:2:2999);

N = 1024;

xk = fft(xn, N);

h = figure;

plot(8001:8000+length(xn), xn);

title('原信号抽取波形');

saveas(h, '原信号抽取波形', 'svg');

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');

%% 思考题 2.

xk = xk(1:2:end);

h = figure;

subplot(2, 1, 1);

plot(linspace(0, 2, N/2), 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/2), angle(xk));

xlabel('\omega (\pi rad/s)');

ylabel('\phi(\omega)');

title('相频响应曲线 \phi(\omega)');

saveas(h, '幅频抽取-频率响应曲线', 'svg');

h = figure;

xn = ifft(xk, N/2);

L = length(xn);

plot(linspace(8000, 9000, L), xn);

xlabel('t');

title('对幅频响应曲线进行抽取后还原的原信号');

saveas(h, '频率抽取-原信号', 'svg');

思考题