DFT 在信号频谱分析中的应用

DFT 在信号频谱分析中的应用

实验目的

  1. 熟悉 DFT 的性质。

    DFT是离散傅里叶变换的缩写，是一种将时域信号转换为频域信号的数学工具。下面是DFT的一些基本性质：

    1. 线性性：DFT是线性的，即它满足叠加原理。如果x1(n)和x2(n)是两个长度为N的离散时间信号，那么它们的DFT可以表示为：

      X(k) = DFT{x(n)}，其中n=0,1,...,N-1
      Y(k) = DFT{x2(n)}，其中n=0,1,...,N-1
      则a1X(k) + a2Y(k)的DFT为：
      Z(k) = a1X(k) + a2Y(k)

    2. 周期性：DFT是周期性的。如果x(n)是一个长度为N的周期信号，即x(n) = x(n+mN)，其中m是任意整数，那么它的DFT可以表示为：

      X(k) = DFT{x(n)}，其中n=0,1,...,N-1
      则X(k)也具有周期性，即X(k) = X(k+mN)

    3. 对称性：如果x(n)是一个实数信号，那么它的DFT具有对称性。具体来说，如果X(k)表示x(n)的DFT，那么有：

      X(k) = X(N-k)，其中表示共轭复数

    下面是用Matlab代码演示DFT这些性质的例子：

    % 线性性举例
x1 = [1 2 3 4];
x2 = [5 6 7 8];
a1 = 2;
a2 = 3;
X1 = fft(x1);
X2 = fft(x2);
Z = a1*X1 + a2*X2;
disp('DFT的线性性：')
disp(Z)

% 周期性举例
x = [1 2 3 4 1 2 3 4];
X = fft(x);
disp('DFT的周期性：')
disp(X)

% 对称性举例
x = [1 2 3 4];
X = fft(x);
disp('DFT的对称性：')
disp(X)
disp(conj(fliplr(X)))

    运行结果为：

    DFT的线性性：
  列 1 至 2

  98.0000 + 0.0000i -10.0000 +10.0000i

  列 3 至 4

 -10.0000 + 0.0000i -10.0000 -10.0000i
DFT的周期性：
  列 1 至 2

  20.0000 + 0.0000i   0.0000 + 0.0000i

  列 3 至 4

  -4.0000 + 4.0000i   0.0000 + 0.0000i

  列 5 至 6

  -4.0000 + 0.0000i   0.0000 + 0.0000i

  列 7 至 8

  -4.0000 - 4.0000i   0.0000 + 0.0000i
DFT的对称性：
  列 1 至 2

  10.0000 + 0.0000i  -2.0000 + 2.0000i

  列 3 至 4

  -2.0000 + 0.0000i  -2.0000 - 2.0000i
  列 1 至 2

  -2.0000 + 2.0000i  -2.0000 + 0.0000i

  列 3 至 4

  -2.0000 - 2.0000i  10.0000 + 0.0000i

    以上代码中，fft函数是Matlab中用于计算DFT的函数。

  2. 深理解信号频谱的概念及性质。

  3. 了解高密度谱与高分辨率频谱的区别。

实验任务与要求

  1. 学习用 DFT 和补零 DFT 的方法来计算信号的频谱。

  2. 用 MATLAB 语言编程来实现，在做课程设计前，必须充分预习课本 DTFT、 DFT 及补零 DFT 的有关概念，熟悉 MATLAB 语言，独立编写程序。

实验仪器设备

MATLAB2022a

实验内容

题目一

用 MATLAB 语言编写计算序列$ x(n)$的 N 点 DFT 的 m 函数文件dft.m。并与 MATLAB 中的内部函数文件fft.m作比较。

题目二

对离散确定信号 x(n)=cos⁡(0.48πn)+cos⁡(0.52πn) 作如下谱分析：

  1. 截取 x(n) 使 x(n) 成为有限长序列 N(0≤n≤N-1), (长度 N 自己选) 写程序计算出 x(n) 的 N 点 DFT X(k), 画出时域序列图 xn∼n 和相应的幅频图 |X(k)|∼k 。
  2. 将 (1) 中 \(x(n)\)补零加长至 M 点, 长度 M 自己选, (为了比较补零长短的影响, M 可以取两次值, 一次取较小的整数, 一次取较大的整数), 编写 程序计算\(x(n)\)的 M 点 DFT, 画出时域序列图和两次补零后相应的 DFT 幅 频图。
  3. 利用补零 DFT 计算 (1) 中 N 点有限长序列 \(x(n)\) 频谱\(X(e^{j\omega})\)并画出相应 的幅频图 \(|X(e^{j\omega})|\sim\omega\)。

题目三

  1. 研究高密度谱与高分辨率频谱。对连续确定信号 \(x_a(t)=\cos(2\pi\cdot6.5\times10^3t)+\cos(2\pi\cdot7\times10^3t)+\cos(2\pi\cdot9\times10^3t)\)，以采样频率 \(f_s=32\text{kHz}\) 对信号 \(x(t)\) 采样得离散信号 \(x(n)\)，分析下列三种情况的幅频特性：

    1. 采集数据 \(x(n)\) 长度取 \(N=16\) 点，编写程序计算出 \(x(n)\) 的 \(16\) 点 DFT \(X(k)\)，并画出相应的幅频图 \(|X(k)|\sim k\)。
    2. 采集数据 \(x(n)\) 长度 \(N=16\) 点，补零加长至 \(M\) 点（长度 \(M\) 自己选），利用补零 DFT 计算 \(x(n)\) 的频谱 \(X_1(e^{j\omega})\) 并画出相应的幅频图 \(|X_1(e^{j\omega})|\sim\omega\)。
    3. 采集数据 \(x(n)\) 长度取为 \(M\) 点（注意不是补零至 \(M\)），编写程序计算出 \(M\) 点采集数据 \(X_2(e^{j\omega})\) 的频谱 \(X_2(e^{j\omega})\) 并画出相应的幅频图 \(|X_2(e^{j\omega})|\sim\omega\)。

实验设计

题目一

以下是用MATLAB语言编写计算序列x(n)的N点DFT的m函数文件dft.m：

function X = dft(x, N)
% Computes the N-point DFT of the sequence x
% Inputs:
%   x: input sequence
%   N: DFT length
% Output:
%   X: DFT sequence

n = 0:N-1; % Time indices
k = 0:N-1; % Frequency indices
WN = exp(-1j*2*pi/N); % Twiddle factor

nk = n'*k; % Outer product of n and k
WNnk = WN .^ nk; % Twiddle factor matrix

X = x * WNnk; % Compute DFT

下面是一个使用dft.m函数的示例：

N = 8; % DFT length
x = [1 2 3 4 4 3 2 1]; % Input sequence

% Compute DFT using dft.m
X = dft(x, N);

% Compute DFT using built-in fft function
X_fft = fft(x, N);

% Compare results
disp(norm(X - X_fft));

上述代码中，我们首先定义了一个函数dft，用于计算输入序列x的N点DFT。该函数使用了矩阵乘法的方法来计算DFT。具体来说，我们首先生成时间和频率指数向量n和k，然后计算n和k的外积，得到一个大小为NxN的矩阵，其中第(i,j)个元素为WN^(i*j)，其中WN是旋转因子。最后，我们将输入序列x乘以这个矩阵，得到DFT序列X。

我们还编写了一个示例程序，使用dft.m函数和MATLAB内置的fft函数来计算DFT，并比较它们的结果。在这个示例中，我们使用了一个长度为8的输入序列，但可以使用任何长度的序列和任何DFT长度来调用dft函数。

下面是示例程序的输出：

2.2414e-14

这表明dft函数和fft函数计算的DFT结果非常接近，差异非常小。因此，我们可以得出结论，dft函数可以与MATLAB内置的fft函数相媲美。

题目二

以下是 MATLAB 代码和结果：

  1. 截取 x(n) 使 x(n) 成为有限长序列 N，计算 N 点 DFT：
% Part 1: calculate N-point DFT of x(n)

N = 64; % choose N as 64
n = 0:N-1;
x = cos(0.48*pi*n) + cos(0.52*pi*n);
X = fft(x, N);

% plot time-domain sequence
subplot(2,1,1);
stem(n, x);
xlabel('n');
ylabel('x(n)');
title('time-domain sequence');

% plot magnitude spectrum
subplot(2,1,2);
stem(linspace(0, 1, N), abs(X));
xlabel('k/N');
ylabel('|X(k)|');
title('magnitude spectrum, N=64');

运行程序后可以得到时域序列图和相应的幅频图如下所示：

img

可以看到，时域序列图中包含了两个正弦波分量，频率分别为 0.48π 和 0.52π。在幅频图中，可以看到有两个峰对应于这两个正弦波分量。

  1. 将 x(n) 补零加长至 M 点，计算 M 点 DFT，比较不同补零长度的效果：
% Part 2: calculate M-point DFT of x(n) with zero-padding

M1 = 256; % choose M1 as 256
M2 = 1024; % choose M2 as 1024
xn1 = [x zeros(1, M1-N)];
xn2 = [x zeros(1, M2-N)];
X1 = fft(xn1, M1);
X2 = fft(xn2, M2);

% plot time-domain sequence
figure;
subplot(2,1,1);
stem(0:M1-1, xn1);
xlabel('n');
ylabel('x1(n)');
title('time-domain sequence, M=256');

subplot(2,1,2);
stem(0:M2-1, xn2);
xlabel('n');
ylabel('x2(n)');
title('time-domain sequence, M=1024');

% plot magnitude spectrum with different zero-padding lengths
figure;
subplot(2,1,1);
stem(linspace(0, 1, M1), abs(X1));
xlabel('k/M');
ylabel('|X1(k)|');
title('magnitude spectrum with zero-padding, M=256');

subplot(2,1,2);
stem(linspace(0, 1, M2), abs(X2));
xlabel('k/M');
ylabel('|X2(k)|');
title('magnitude spectrum with zero-padding, M=1024');

运行程序后可以得到补零前后的时域序列图和幅频图如下所示：

img

可以看到，对于 M=256 的情况，补零后的幅频图与原始幅频图形状基本相同，但分辨率更高，零值点更多；而对于 M=1024 的情况，补零后的幅频图中出现了更多的零值点，且峰更加明显，但是峰之间的分辨率变得更低。

  1. 利用补零 DFT 计算 N 点有限长序列 x(n) 频谱 X(e^jω )：
% Part 3: calculate frequency spectrum of x(n) using zero-padding DFT

w = linspace(0, 2*pi, 1000);
Xw = fft(xn1, M1);
Xew = fft(xn1, M2);

% plot magnitude spectrum
figure;
subplot(2,1,1);
plot(w, abs(interp1(linspace(0, 2*pi, M1), Xw, w, 'spline')));
xlabel('w');
ylabel('|X(e^jw)|');
title('magnitude spectrum with zero-padding, M=256');

subplot(2,1,2);
plot(w, abs(interp1(linspace(0, 2*pi, M2), Xew, w, 'spline')));
xlabel('w');
ylabel('|X(e^jw)|');
title('magnitude spectrum with zero-padding, M=1024');

运行程序后可以得到幅频图如下所示：

img

img

可以看到，频谱中有两个峰对应于信号中的两个正弦波分量，与 DFT 时域序列长度为 N 时的结果相同。但是，频谱分辨率更高，零值点更多，这是由于补零操作增加了数据点的数量。

题目三

以下是实现要求的 MATLAB 代码：

% define the continuous signal x_a(t)
t = 0:1/32000:0.02; % define the time vector with a sampling frequency of 32 kHz
xa = cos(2*pi*6500*t) + cos(2*pi*7000*t) + cos(2*pi*9000*t);

% sample the continuous signal x_a(t) to obtain the discrete signal x(n)
x = xa(1:16);

% calculate the 16-point DFT of x(n) and plot the magnitude spectrum
X = fft(x);
figure;
stem(abs(X));
xlabel('k');
ylabel('|X(k)|');
title('Magnitude Spectrum of 16-point DFT');

% zero-pad x(n) to obtain a signal of length M
M = 128;
xzp = [x zeros(1,M-16)];

% calculate the DFT of xzp(n) and plot the magnitude spectrum
X1 = fft(xzp);
omega = 2*pi/M*(0:M-1);
figure;
plot(omega, abs(X1));
xlabel('\omega');
ylabel('|X_1(e^{j\omega})|');
title('Magnitude Spectrum of Zero-Padded DFT');

% sample the continuous signal x_a(t) to obtain a signal of length M
xaM = cos(2*pi*6500*(0:1/32000:(M-1)*1/32000)) + cos(2*pi*7000*(0:1/32000:(M-1)*1/32000)) + cos(2*pi*9000*(0:1/32000:(M-1)*1/32000));

% calculate the DFT of xaM(n) and plot the magnitude spectrum
X2 = fft(xaM);
omega = 2*pi/M*(0:M-1);
figure;
plot(omega, abs(X2));
xlabel('\omega');
ylabel('|X_2(e^{j\omega})|');
title('Magnitude Spectrum of Sampled DFT');

运行程序后可以得到三张幅频特性图：

  1. 采集数据 \(x(n)\) 长度取 \(N=16\) 点的幅频特性图：
    img

  2. 采集数据 \(x(n)\) 长度 \(N=16\) 点，补零加长至 \(M=128\) 点的幅频特性图：
    img

  3. 采集数据 \(x(n)\) 长度取为 \(M=128\) 点的幅频特性图：
    img

可以根据需要修改采样频率 \(f_s\)、采集数据长度 \(N\) 和加长后的数据长度 \(M\)。

实验结论探讨及分析

思考题回答

假设实际测得的一段信号的表达式为：

\(x(t)=0.5cos(2\pi f_2t)+0.75cos(2\pi f_2t)\)

其中\(f_1\),\(f_2\)自定

试确定一合适的采样频率 \(f_1\) ，利用 fft 分析该信号的频谱。在信号截短时 要求

1） 使用矩形窗，考虑频谱泄露和频率分辨率等的影响，确定采样数据的长度。
画出信号的时域和频域波形。
2） 使用汉宁窗，确定能够分辨出最小谱峰间距的信号长度，并画出对应的信
号时域和频域波形图。

假设\(f_1=100\)Hz，\(f_2=500\)Hz，我们需要确定合适的采样频率\(f_s\)，使得信号能够被充分采样，并且能够准确地分析出其频谱。

首先，我们需要确定信号的最高频率\(f_{max}\)。根据奈奎斯特采样定理，采样频率应该至少为最高频率的两倍，即\(f_s \geq 2f_{max}\)。在这个例子中，\(f_{max}=500\)Hz，因此我们可以选择\(f_s=1000\)Hz作为采样频率。

接下来，我们将信号进行采样。为了确定采样数据的长度，我们可以考虑频谱泄露和频率分辨率等因素。我们将采用两种窗函数：矩形窗和汉宁窗。

对于矩形窗，我们需要确定采样数据的长度\(N\)。根据频率分辨率的公式\(\Delta f=1/T=N/T_s\)，其中\(T\)为采样数据的总时间，\(T_s\)为采样间隔。我们可以选择\(T=1\)秒，\(T_s=1/f_s=1/1000\)秒，因此\(\Delta f=N/T_s=Nf_s\)。为了分辨两个频率分量，我们需要确保它们之间的距离大于\(\Delta f\)，即\(500-100 > Nf_s\)。因此，我们可以选择\(N=400\)。

对于汉宁窗，我们需要确定能够分辨出最小谱峰间距的信号长度。根据汉宁窗的频谱主瓣宽度公式\(\Delta f=2/T\)，其中\(T\)为采样数据的总时间，我们可以选择\(T=1\)秒，因此\(\Delta f=2\)Hz。为了分辨两个频率分量，我们需要确保它们之间的距离大于\(\Delta f\)，即\(500-100 > 2\)Hz\(N\)。因此，我们可以选择\(N=200\)。

下面是Matlab代码实现：

% 定义信号 x(t)
f1 = 100;
f2 = 500;
t = 0:1/1000:1-1/1000;
x = 0.5*cos(2*pi*f1*t) + 0.75*cos(2*pi*f2*t);

% 矩形窗
N1 = 400;
X1 = fft(x(1:N1).*rectwin(N1)', N1);
f1 = (0:N1-1)*1000/N1;
figure;
subplot(2, 1, 1);
plot(t(1:N1), x(1:N1));
xlabel('Time (s)');
ylabel('Amplitude');
title('Rectangular Window');
subplot(2, 1, 2);
stem(f1, abs(X1));
xlabel('Frequency (Hz)');
ylabel('Magnitude');

% 汉宁窗
N2 = 200;
X2 = fft(x(1:N2).*hann(N2)', N2);
f2 = (0:N2-1)*1000/N2;
figure;
subplot(2, 1, 1);
plot(t(1:N2), x(1:N2));
xlabel('Time (s)');
ylabel('Amplitude');
title('Hanning Window');
subplot(2, 1, 2);
stem(f2, abs(X2));
xlabel('Frequency (Hz)');
ylabel('Magnitude');

代码首先定义了信号\(x(t)\)，然后使用矩形窗和汉宁窗分别对信号进行采样和频谱分析，并绘制了时域和频域波形。可以看出，使用矩形窗时，频谱主瓣宽度较宽，谱峰之间存在重叠，无法准确地分辨出两个频率分量。而使用汉宁窗时，频谱主瓣宽度较窄，谱峰之间无重叠，能够准确地分辨出两个频率分量。

因此，根据我们的分析，使用汉宁窗时，能够分辨出最小谱峰间距的信号长度为\(N=200\)，对应的频谱图显示出两个明显的峰值。****

posted on 2023-09-04 19:24  白茶木  阅读(222)  评论(0)  编辑  收藏  举报
Powered by:
博客园
Copyright © 2024 白茶木
Powered by .NET 9.0 on Kubernetes