信号处理实验——分析采样间隔和采样长度对FFT的影响

最新推荐文章于 2024-09-20 09:03:12 发布

weixin_45951391

最新推荐文章于 2024-09-20 09:03:12 发布

阅读量1.9k

点赞数 28

文章标签：信号处理人工智能

本文链接：https://blog.csdn.net/weixin_45951391/article/details/136484836

版权

实验目的

分析采样间隔和采样长度对FFT 的影响
分析栅栏效应对FFT的影响

实验原理

由于非周期信号（周期信号截短后也为非周期信号）具有连续的频谱，而用FFT只能计算出其离散频谱，即连续频谱中的若干点，这就好像通过栅栏的缝隙观看另一边，只能在离散点处看到真实的频谱，这种现象称为栅栏效应。此时，FFT像一个“栅栏”用来观察周期信号的离散频谱，而这些离散谱恰恰是处于透过栅栏观察不到的部分

实验方法与内容

已知余弦信号如下 $x(t)=cos(2\pi Ft)$ 对该信号进行傅里叶变换，用FFT程序分析正弦结果，分别在以下情况下进行，并且分析比较结果。

编程运行：

仿真代码：

% 定义参数
F = 50; % 信号频率
Ns = [32, 32, 32, 32, 64]; % 采样长度数组
Ts = [0.000625, 0.005, 0.0046875, 0.004, 0.000625]; % 采样间隔数组

% 循环绘制FFT分析结果
for i = 1:length(Ns)
    [Xk, m] = fft_analysis(F, Ns(i), Ts(i));
    subplot(5, 1, i);
    stem(m, abs(Xk));
    title(['F=', num2str(F), ', N=', num2str(Ns(i)), ', T=', num2str(Ts(i))]);
    xlabel('Frequency (Hz)');
    ylabel('Magnitude');
end

% 自定义函数：FFT分析
function [Xk, m] = fft_analysis(F, N, T)
    fs = 1/T; % 采样率
    n = 0:N-1;
    xn = cos(2*pi*F*n/fs); % 信号模型
    Xk = fft(xn); % 快速傅里叶变换
    m = 0:N-1; % 频率序列
end

运行结果：

虽然这五种情形下均满足抽样定理，但FFT的结果却并不如预期所愿为两束脉冲，必须对其进行分析。根据余弦函数的DFT通式 $cos(\frac{2\pi l}{N}n )\Leftrightarrow N[\delta (k-l)+\delta (k+l)]$ 可以计算 $cos(2\pi FTn )$ 在这5种不同情况下的DFT，可以分别得到 $l=FTn$ 的值，即冲激所在的位置：

由上图可以看出，对于1、2、5三种情况， $\delta$ 为整数，与实验结果图完美契合。

但对于3、4这两种情况， $\delta$ 为小数，而FFT的 k值只能在整数点取值，猜测是受栅栏效应的影响， $\delta$ 点的响应值无法看见。根据此猜想，下面以第2种情况和第4种情况为例，即T=0.005和T=0.004，对其进行进一步分析。

运行代码：

n = 0:31;
k = 0:999; % 调整以匹配FFT长度
w = (pi/500)*k;

% 自定义FFT分析函数（假设在脚本的其他地方已定义）
function [Xk, m] = fft_analysis(F, N, T)
    fs = 1/T; % 采样率
    n = 0:N-1;
    xn = cos(2*pi*F*n/fs); % 信号模型
    Xk = fft(xn); % 快速傅立叶变换
    m = 0:N-1; % 频率区间
end

% F=50, N=32, T=0.005
s1 = cos(2*pi*50*0.005*n);
DFT_matrix = exp(-1i*pi/500).^(n' * k); % 直接DFT计算矩阵
S1 = s1 * DFT_matrix;
[Xk2, m2] = fft_analysis(50, 32, 0.005);

subplot(2,1,1)
plot(w/pi*16, abs(S1)); hold on; % 修正缩放以便比较
stem(m2, abs(Xk2), 'r'); % FFT用不同颜色
title('F=50, N=32, T=0.005'); xlabel('频率 (π/16单位)'); ylabel('幅度');
axis([0 32 0 20]); hold off;

% F=50, N=32, T=0.004
s2 = cos(2*pi*50*0.004*n);
S2 = s2 * DFT_matrix; % 使用更新的信号重新使用DFT矩阵
[Xk4, m4] = fft_analysis(50, 32, 0.004);

subplot(2,1,2)
plot(w/pi*16, abs(S2)); hold on; % 一致缩放
stem(m4, abs(Xk4), 'r'); % FFT用不同颜色
title('F=50, N=32, T=0.004'); xlabel('频率 (π/16单位)'); ylabel('幅度');
axis([0 32 0 20]); hold off;

运行结果：

作出这两种情况下的DTFT，并分别与原FFT的图像进行比较。对比结果图如上所示。很明显，受栅栏效应的影响，FFT只是DTFT在2π/N时的抽样点。对于T=0.005，其FFT中纵坐标为0的众多点，只是恰好取到了其DTFT的零点，而FFT中独立的两个脉冲只是恰好取到了DTFT的最大值。同样的，由于T取值的特殊性，对于T=0.004，其FFT并不能恰好取到DTFT的零点，也就出现了FFT中起伏波动的频率点。（也未必刚好可以抽样到峰值，比如可以将T改为0.0039）

为了改善栅栏效应造成的影响，尝试在序列后补零。

以下分别为补1倍长度的0和3倍长度的0的源代码和对比结果。

1）运行代码：

n = 0:31; % 信号序列的时间索引
k = 0:999; % 频率分点，用于DFT计算
w = (pi/500) * k; % 频率变量，用于绘制频谱

% 序列后增补了相等长度的0，意味着总长度变为原来的2倍
% F=50, N=32, T=0.005
subplot(2,1,1);
s1 = cos(2*pi*50*0.005*n); % 生成信号
DFT_matrix = exp(-1i*pi/500).^(n' * k); % 计算DFT矩阵
S1 = s1 * DFT_matrix; % 应用DFT
[Xk2, m2] = fft_analysis(50, 32, 0.005); % 使用自定义FFT分析函数
plot(w/pi*16, abs(S1)); % 绘制连续DFT结果
axis([0 32 0 20]); hold on;
stem(m2/2, abs(Xk2), 'r'); % 绘制FFT结果，红色以区分
title('F=50, N=32, T=0.005');
xlabel('频率 (单位：π/16)');
ylabel('幅度');
hold off;

% F=50, N=32, T=0.004
subplot(2,1,2);
s2 = cos(2*pi*50*0.0039*n); % 第二个信号模型，注意T的改变
S2 = s2 * DFT_matrix; % 重用DFT矩阵
[Xk4, m4] = fft_analysis(50, 32, 0.0039); % 自定义FFT分析
plot(w/pi*16, abs(S2)); % 绘制连续DFT结果
axis([0 32 0 20]); hold on;
stem(m4/2, abs(Xk4), 'r'); % 绘制FFT结果，红色以区分
title('F=50, N=32, T=0.0039');
xlabel('频率 (单位：π/16)');
ylabel('幅度');
hold off;

function [Xk, m] = fft_analysis(F, N, T)
    fs = 1/T; % 计算采样率
    n = 0:N-1; % 时间索引
    xn = [cos(2*pi*F*n/fs), zeros(1, N)]; % 生成信号并增补0
    Xk = fft(xn); % 应用快速傅里叶变换
    m = 0:(2*N-1); % 频率索引范围扩大为2倍
end

仿真结果：

2）运行代码：

n = 0:31; % 时间序列索引
k = 0:999; % 频率点索引，用于计算DFT
w = (pi/500) * k; % 频率变量，用于绘图

% 对序列进行3倍长度的0填充
% F=50, N=32, T=0.005
subplot(2,1,1);
s1 = cos(2*pi*50*0.005*n); % 生成信号
DFT矩阵 = exp(-1i*pi/500).^(n' * k); % 计算DFT矩阵
S1 = s1 * DFT矩阵; % 应用DFT
[Xk2, m2] = fft_analysis(50, 32, 0.005); % 使用自定义FFT分析函数
plot(w/pi*16, abs(S1)); % 绘制DFT结果
axis([0 32 0 20]); hold on;
stem(m2/4, abs(Xk2), 'r'); % 绘制FFT结果，红色以区分
title('F=50, N=32, T=0.005');
xlabel('频率 (单位：π/16)');
ylabel('幅度');
hold off;

% F=50, N=32, T=0.0039
subplot(2,1,2);
s2 = cos(2*pi*50*0.0039*n); % 第二个信号模型，注意T的改变
S2 = s2 * DFT矩阵; % 重用DFT矩阵
[Xk4, m4] = fft_analysis(50, 32, 0.0039); % 自定义FFT分析
plot(w/pi*16, abs(S2)); % 绘制DFT结果
axis([0 32 0 20]); hold on;
stem(m4/4, abs(Xk4), 'r'); % 绘制FFT结果，红色以区分
title('F=50, N=32, T=0.0039');
xlabel('频率 (单位：π/16)');
ylabel('幅度');
hold off;

% 自定义FFT分析函数
function [Xk, m] = fft_analysis(F, N, T)
    fs = 1/T; % 采样率
    n = 0:N-1; % 时间序列索引
    xn = [cos(2*pi*F*n/fs), zeros(1, 3*N)]; % 生成信号并进行3倍长度的0填充
    Xk = fft(xn); % 执行FFT
    m = 0:(4*N-1); % 调整频率点索引范围，考虑到了0填充
end

仿真结果：