频域分析实践介绍

2023-05-22 15:42| 来源: 网络整理| 查看: 265

FFT 的幅值和相位信息

信号的频域表示携带每个频率下信号幅值和相位的有关信息。因此，FFT 计算的输出是复数。复数有实部和虚部，满足。的幅值计算为，的相位计算为。使用 MATLAB 函数 abs 和 angle 可以分别得到任意复数的幅值和相位。

使用一个音频示例来了解信号的幅值和相位携带哪些信息。为此，加载一个包含 15 秒原声吉他音乐的音频文件。音频信号的采样率为 44.1 kHz。

Fs = 44100; y = audioread('guitartune.wav');

使用 fft 观察信号的频率成分。

NFFT = length(y); Y = fft(y,NFFT); F = ((0:1/NFFT:1-1/NFFT)*Fs).';

FFT 的输出是一个复数向量，其中包含关于信号频率成分的信息。幅值显示频率分量相对于其他分量的强度。相位显示所有频率分量按时间对齐的情况。

绘制信号频谱的幅值和相位分量。幅值以对数刻度 (dB) 方便地绘制。使用 unwrap 函数展开相位，以便我们可以看到频率的连续函数。

magnitudeY = abs(Y); % Magnitude of the FFT phaseY = unwrap(angle(Y)); % Phase of the FFT helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)

您可以对频域向量 Y 应用傅里叶逆变换以还原时间信号。'symmetric' 标志告知 ifft 您正在处理实数值时间信号，因此它会将由于计算中的数值不准确而在逆变换上出现的较小虚部归零。请注意，原始时间信号 y 和还原后的信号 y1 实际上是相同的（其差值的范数的数量级为 1e-14）。两者之间的微小差值也是由上面提到的数值不准确造成的。播放并聆听未变换的信号 y1。

y1 = ifft(Y,NFFT,'symmetric'); norm(y-y1) ans = 3.7851e-14 hplayer = audioplayer(y1, Fs); play(hplayer);

要查看更改信号幅值响应的效果，请直接从 FFT 输出中去除 1 kHz 以上的频率分量（通过使幅值等于零），并通过聆听来了解此操作对音频文件声音的影响。去除信号的高频分量称为低通滤波。

Ylp = Y; Ylp(F>=1000 & F N 可有效地对频率点进行插值，从而呈现更详细的频谱图（这是通过在时间信号的末尾追加 NFFT-N 个零点并选取零填充向量的 NFFT 点 FFT 来实现的）。

开环和闭环加速度功率谱表明，当控制系统处于活动状态时，加速度功率谱在 4 和 11 dB 之间减小。最大衰减发生在大约 23.44 kHz 处。11 dB 的降幅意味着振动功率按因子 12.6 降低。总功率从 0.1670 瓦特降低到 0.059 瓦特，因子为 2.83。

load quakevibration.mat Fs = 1e3; % sample rate NFFT = 512; % number of FFT points segmentLength = 64; % segment length % open loop acceleration power spectrum [P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power'); % closed loop acceleration power spectrum P1_CL = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power'); helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),... 'Frequency in Hz','Acceleration Power Spectrum in dB',... 'Resolution bandwidth = 15.625 Hz',{'Open loop', 'Closed loop'},[0 100])

您正在分析振动数据并知道振动有循环行为。那么，为什么上面显示的频谱图不包含任何典型循环行为的尖锐频谱线？不显示这些线，是否因为它们无法用包含 64 个点的段长度获得的分辨率来解析？提高频率分辨率，看看是否存在以前无法解析的频谱线。为此，请将 pwelch 函数中使用的数据段长度增加到 512 个点。这会产生 Fs/512 = 1.9531 Hz 的新分辨率。在本例中，FFT 平均值的数目减少到 floor(10e3/512) = 19 个。显然，使用 pwelch 时，平均值的数目和频率分辨率之间存在权衡。保持 FFT 点数等于 512。

NFFT = 512; % number of FFT points segmentLength = 512; % segment length [P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power'); P1_CL = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power'); helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),... 'Frequency in Hz','Acceleration Power Spectrum in dB',... 'Resolution bandwidth = 1.95 Hz',{'Open loop', 'Closed loop'},[0 100])

请注意，频率分辨率的提高使您能够观测到开环频谱上的三个峰值和闭环频谱上的两个峰值。这些峰值以前是无法解析的。开环频谱上的峰值之间的间隔约为 11 Hz，小于用长度为 64 的段获得的频率分辨率，但大于用长度为 512 的段获得的分辨率。振动的循环行为现在可见。主振动频率为 5.86 Hz，等间距的频率峰值表明它们是在谐波上相关。虽然已观测到控制系统可降低振动的总功率，但是，更高分辨率频谱表明控制系统的另一个效果是使谐波分量在 17.58 Hz 处出现陷波。因此，控制系统不仅减少振动，还使其更接近正弦波。

请注意，频率分辨率由信号点的数量确定，而不是由 FFT 点的数量确定。增加 FFT 点的数量会对频率数据进行插值，以提供更多频谱细节，但不会提高分辨率。

总结

在此示例中，您学习了如何使用 fft、ifft、periodogram、pwelch 和 bandpower 函数对信号执行频域分析。您已经了解 FFT 的复/实性质，以及频谱的幅值和相位中包含哪些信息。您看到了在分析信号的周期性时使用频域数据的优势。您已经掌握如何计算总功率或含噪信号的特定频带的功率。您已经了解如何通过提高频谱的频率分辨率来观测近距频率分量，并了解频率分辨率和频谱求平均值之间的权衡。