主要内容

频域分析实用导论

此示例演示如何执行和解释基本频域信号分析。该示例讨论了使用信号频域表示法与时域表示法的优势,并使用模拟数据和真实数据说明了基本概念。该示例回答了一些基本问题,如:幅值和pha的含义是什么FFT的se?我的信号是周期性的吗?我如何测量功率?这个频带中有一个或多个信号吗?

频域分析是信号处理应用中最重要的工具。频域分析广泛应用于通信、地质、遥感和图像处理等领域。时域分析显示信号如何随时间变化,而频域分析显示信号能量如何分布在频率范围内。频域表示还包括必须应用于每个频率分量的相移信息,以便通过所有单独频率分量的组合恢复原始时间信号。

信号可以用一对数学运算符在时域和频域之间进行转换,称为变换。一个例子是傅里叶变换,它将一个函数分解成一个(可能是无限的)正弦波频率分量的和。频率分量的“频谱”是信号的频域表示。傅里叶反变换将频域函数转换回时间函数。的fftifftMATLAB中的函数允许您分别计算信号的离散傅里叶变换(DFT)和这个变换的逆。

FFT的幅度和相位信息

信号的频域表示法携带有关信号在每个频率上的幅度和相位的信息。这就是为什么FFT计算的输出是复杂的。一个复数,x美元,有实部,x_r美元,虚部,x_i美元,这样$x = x_r + jx_i$.的大小x美元是计算$ \√6 {(x_r ^ 2 + x_i ^ 2)} $,以及x美元是计算$\arctan{(x_i/x_r)}$.你可以使用MATLAB函数防抱死制动系统分别求任意复数的幅值和相位。

以音频为例,了解信号的幅度和相位所承载的信息。为此,加载一个包含15秒原声吉他音乐的音频文件。音频信号的采样率为44.1 kHz。

Fs = 44100;y = audioread (“guitartune.wav”);

使用fft观察信号的频率含量。

NFFT=长度(y);y=fft(y,NFFT);F=((0:1/NFFT:1-1/NFFT)*Fs)。”;

FFT的输出是一个包含信号频率内容的信息的复向量。振幅告诉你相对于其他分量频率分量的强度。相位告诉你所有的频率分量是如何在时间上对齐的。

绘制信号频谱的幅度和相位分量。星等可以方便地用对数标度(dB)表示。阶段使用打开函数,这样我们就能看到频率的连续函数。

magnitudeY = abs (Y);% FFT的幅度相位=展开(角度(Y));%相位的FFThelperFrequencyAnalysisPlot1 (F, magnitudeY phaseY NFFT)

你可以对频域向量Y进行傅里叶反变换,来恢复时间信号。'symmetric'标志告诉我们ifft你处理的是一个实值时间信号所以它会使逆变换中出现的虚分量归零因为计算中的数值不准确。注意,原始的时间信号y和恢复的信号y1实际上是相同的(它们的差的范数在1e-14的数量级上)。两者之间的微小差异也是由于上面提到的数值误差造成的。播放并收听未转换的信号y1。

日元=传输线(Y, NFFT,“对称”);规范(y-y1)
ans = 3.7851 e-14
hplayer = audioplayer(y1, Fs);玩(hplayer);

为了看到改变信号的幅度响应的影响,直接从FFT输出中去除1 kHz以上的频率成分(通过使幅度等于零),并听这对音频文件的声音的影响。去除信号中的高频成分称为低通滤波。

Ylp = Y;Ylp(F>=1000 & F<=Fs-1000) = 0;helperFrequencyAnalysisPlot1 (F、abs (Ylp),打开(角(Ylp)), NFFT,...“1 kHz以上的频率分量已归零”

使用以下方法将滤波后的信号恢复到时域ifft

ylp =传输线(ylp,“对称”);

播放信号。你仍然可以听到旋律,但听起来就像你捂住耳朵一样(你这样做时会过滤高频声音)。即使吉他产生的音符在400到1千赫之间,当你在琴弦上弹奏音符时,琴弦也会以基频的倍数振动。这些被称为谐波的高频成分赋予吉他独特的音色。当你移除它们时,你会使声音看起来“不透明”。

hplayer = audioplayer(ylp, Fs);玩(hplayer);

信号的相位有关于歌曲的音符何时出现的重要信息。为了说明相位对音频信号的重要性,可以通过取每个频率分量的幅值来完全去除相位信息。注意,这样做可以保持幅度响应不变。

%取信号的每个FFT分量的幅值Yzp = abs (Y);helperFrequencyAnalysisPlot1 (F、abs (Yzp),打开(角(Yzp)), NFFT, [],...“相位已设置为零”

把信号放回时域,然后播放音频。你根本认不出原来的声音。幅度响应是相同的,这次没有去掉频率成分,但是音符的顺序完全消失了。现在信号由一组在时间等于零时对齐的正弦波组成。一般来说,滤波引起的相位畸变会使信号损坏到无法识别的地步。

yzp=ifft(yzp,“对称”);hplayer = audioplayer(yzp, Fs); / /播放玩(hplayer);

发现信号周期的研究

信号的频域表示允许你观察信号的几个特征,这些特征要么是不容易看到的,要么是当你在时域观察信号时根本看不见的。例如,当您寻找信号的循环行为时,频域分析就变得非常有用。

某办公楼温度循环行为分析

考虑一组办公楼在冬季的温度测量。在大约16.5周的时间里,每30分钟测量一次。看看时域数据,时间轴按周缩放。这些数据是否存在周期性行为?

负载officetemp.matFs = 1 / (60 * 30);%采样率为每30分钟1个样本t =(0:长度(临时)1)/ Fs;helperFrequencyAnalysisPlot2 (t /(60 * 60 * 24 * 7),温度,...“在周”“温度(华氏)”

通过观察时域信号,几乎不可能知道办公室温度是否存在任何循环行为。然而,温度的循环行为变得明显,如果我们看它的频域表示。

获得信号的频域表示。如果你将FFT输出的幅值用频率轴按周期/周缩放,你可以看到有两条谱线,显然比任何其他频率分量都要大。一条谱线为1周期/周,另一条谱线为7周期/周。这是有意义的,因为数据来自一个7天日历上的温度控制建筑。第一条谱线表明,建筑温度遵循每周循环,周末温度较低,工作日温度较高。第二行表示也有一个昼夜循环,夜间温度较低,白天温度较高。

NFFT =长度(临时);% FFT点数F = (0: 1/NFFT: 1/2-1/NFFT)*F;%频率矢量TEMP = fft(临时、NFFT);临时(1)= 0;%删除直流组件以更好地显示helperFrequencyAnalysisPlot2 (F * 60 * 60 * 24 * 7, abs(临时(1:NFFT / 2)),...‘频率(周期/周)’“震级”10 [] [], [0])

测量能力

周期图函数计算信号的FFT并对输出进行归一化,以获得功率谱密度(PSD)或功率谱,从中可以测量功率。PSD描述了时间信号的功率如何随频率分布,它的单位是瓦/赫兹。通过将PSD的每个点在该点定义的频率区间上积分(即在PSD的分辨率带宽上),可以计算功率谱。功率谱的单位是瓦。您可以直接从功率谱读取功率值,而不必在一个区间内积分。注意,PSD和功率谱是真实的,所以它们不包含任何相位信息。

非线性功率放大器输出谐波的测量

加载在具有三阶畸变形式的功率放大器的输出处测量的数据$v_o = v_i + 0.75 v_i^2 + 0.5 v_i^3$,在那里v_o美元是输出电压和v_i美元为输入电压。采集的数据采样率为3.6 kHz。输入v_i美元由一个60赫兹、振幅统一的正弦信号组成。由于非线性失真的性质,你应该期望放大器输出信号包含一个直流分量,一个60hz分量,和第二次和第三次谐波在120和180赫兹。

负载放大器输出的3600个样本,计算功率谱,并将结果绘制成对数尺度(分贝-瓦或dBW)。

负载ampoutput1.matFs = 3600;NFFT =长度(y);当你通过“功率”标志输入时计算功率谱[P F] =周期图(y, [], NFFT Fs,“权力”);helperFrequencyAnalysisPlot2 (F, 10 * log10 (P),“频率单位为Hz”...的功率谱(瓦分贝),[],[],[-0.5 200])

功率谱图显示了四个预期峰值中的三个,分别为直流、60和120 Hz。它还显示了多个杂散峰值,这些杂散峰值必须由信号中的噪声引起。请注意,180 Hz谐波完全隐藏在噪声中。

测量可见预期峰值的功率:

PdBW=10*log10(P);在直流时的功率dBW=PdBW(F=0)%瓦分贝[peakPowers_dBW,peakFreqIdx]=FindPeak(PdBW,“minpeakheight”, -11);peakFreqs_Hz = F(peakFreqIdx
power_at_DC_dBW = -7.8873 peakFreqs_Hz = 60 120 peakPowers_dBW = -0.3175 -10.2547

改进噪声信号的功率测量

如上图所示,周期图显示了几个与感兴趣的信号无关的频率峰值。频谱看起来很嘈杂。原因是你只分析了一个短实现的噪声信号。重复这个实验几次,然后取平均值,就可以去除伪谱峰,得到更精确的功率测量结果。可以使用pwelch函数。该函数将取一个大的数据向量,将其分解为指定长度的较小分段,计算分段中尽可能多的周期图,并对它们进行平均。随着可用片段数量的增加,pwelch函数将产生更平滑的功率谱(方差更小),功率值更接近预期值。

加载一个更大的观察,包括500e3点的放大器输出。保持执行fft的点数为3600,使floor(500e3/3600) = 138 fft取平均值,得到功率谱。

负载ampoutput2.matSegmentLength = NFFT;当你通过“功率”标志输入时计算功率谱[P F] = pwelch (y)的(SegmentLength, 1), 0, NFFT, Fs,“权力”);helperFrequencyAnalysisPlot2 (F, 10 * log10 (P),“频率单位为Hz”...的功率谱(瓦分贝),[],[],[-0.5 200])

从情节上看,,pwelch有效地消除了由噪声引起的所有杂散频率峰值。180 Hz处掩埋在噪声中的频谱分量现在可见。平均值消除了频谱的差异,从而有效地产生更精确的功率测量。

测量总平均功率和频带功率

测量时域信号的总平均功率是一项简单而常见的任务。对于放大器输出信号y,在时域中计算总平均功率为:

压水式反应堆= (y ^ 2)之和/长度(y)%的美国瓦茨
压水堆=8.1697

在频域,总平均功率是信号所有频率分量的功率之和。pwr1的值由信号功率谱中所有可用频率分量的和组成。该值与上面使用时域信号计算的pwr值一致:

pwr1 = (P)和%的美国瓦茨
pwr1 = 8.1698

但是,如果您想要测量一个频带上的总功率,该怎么办呢?你可以使用bandpower函数来计算任意频带上的功率。您可以将时域信号直接作为该函数的输入来获得指定频带上的功率。在这种情况下,函数将用周期图法估计功率谱。

计算50hz ~ 70hz范围内的功率。结果将包括60hz的功率加上感兴趣频带上的噪音功率:

pwr_band =带宽功率(y,Fs,[50 70]);pwr_band_dBW = 10 * log10 (pwr_band)%瓦分贝
pwr_band_dBW = 0.0341

如果想控制用于测量频带内功率的功率谱的计算,可以将PSD向量传递给bandpower作用例如,您可以使用pwelch函数来计算PSD,并确保噪声效果的平均:

当您指定“psd”选项时,将计算%功率谱密度(PSD, F) = pwelch (y)的(SegmentLength, 1), 0, NFFT, Fs,“psd”); pwr_波段1=波段功率(PSD,F,[50 70],“psd”);pwr_band_dBW1 = 10 * log10 (pwr_band1)%瓦分贝
pwr_band_dBW1 = 0.0798

发现谱组件

一个信号可能由一个或多个频率成分组成。观察所有光谱成分的能力取决于分析的频率分辨率。功率谱的频率分辨率或分辨率带宽定义为R = Fs/N,其中N为信号观测长度。只有频率大于频率分辨率的频谱成分才会被分辨出来。

某建筑地震振动控制系统分析

主动质量驱动器(AMD)控制系统用于降低建筑物在地震下的振动。一个有源质量驱动器放置在建筑的顶层,根据建筑楼层的位移和加速度测量,一个控制系统向驱动器发送信号,使质量移动以衰减地面扰动。在地震条件下,对一个三层测试结构的一层进行了加速度测量。测量是在没有主动质量驱动控制系统(开环条件)和有主动控制系统(闭环条件)的情况下进行的。

加载加速度数据,计算一楼加速度的功率谱。数据向量的长度为10e3,采样率为1khz。使用pwelch以长度为64个数据点的段获得底部(10e3/64) = 156 FFT平均值和Fs/64 = 15.625 Hz的分辨率带宽。如前所述,平均降低了噪声影响并产生更精确的功率测量。使用512 FFT点。使用NFFT > N有效地插值频率点,绘制更详细的频谱图(这是通过在时间信号的末尾附加NFFT-N零并对填充的零向量取NFFT点FFT来实现的)。

开环和闭环加速度功率谱表明,当控制系统处于活动状态时,加速度功率谱在4到11 dB之间减小。最大衰减发生在23.44 kHz左右。11 dB的减小意味着振动功率减小了12.6倍。总功率从0.1670瓦减小到0.059瓦s、 系数为2.83。

负载quakevibration.matFs = 1 e3;%采样率NFFT = 512;% FFT点数segmentLength = 64;%区段长度%开环加速度功率谱[P1_OL F] = pwelch (gfloor1OL (segmentLength, 1), 0, NFFT, Fs,“权力”);%闭环加速度功率谱P1_CL = pwelch (gfloor1CL (segmentLength, 1), 0, NFFT, Fs,“权力”);helperFrequencyAnalysisPlot2 (F, 10 * log10 (((P1_OL) (P1_CL))),...“频率单位为Hz”“加速度功率谱(dB)”...'分辨率带宽= 15.625 Hz', {“开环”“闭环”}, 100年[0])

您正在分析振动数据,并且您知道振动具有循环行为。那么,上面显示的谱图怎么不包含任何典型循环行为的尖锐谱线呢?也许您缺少这些线是因为它们不能用64点线段长度获得的分辨率进行解析?增加频率分辨率以查看是否存在以前无法解析的谱线。为此,请增加中使用的数据段长度pwelch函数设置为512点。这将产生一个新的分辨率Fs/512=1.9531 Hz。在这种情况下,FFT平均数减少到下限(10e3/512)=19。显然,使用时平均数和频率分辨率之间存在权衡pwelch.保持FFT点数等于512。

NFFT = 512;% FFT点数分段长度=512;%区段长度[P1_OL F] = pwelch (gfloor1OL (segmentLength, 1), 0, NFFT, Fs,“权力”);P1_CL = pwelch (gfloor1CL (segmentLength, 1), 0, NFFT, Fs,“权力”);helperFrequencyAnalysisPlot2 (F, 10 * log10 (((P1_OL) (P1_CL))),...“频率单位为Hz”“加速度功率谱(dB)”...'分辨率带宽= 1.95 Hz', {“开环”“闭环”}, 100年[0])

注意频率分辨率的增加是如何让你在开环频谱上观察到三个峰和在闭环频谱上观察到两个峰的。这些峰以前是无法分解的。开环谱峰间距约为11 Hz,小于64段时的频率分辨率,但大于512段时的频率分辨率。现在可以看到振动的循环特性了。主振动频率为5.86 Hz,等频峰表明两者具有谐波关系。虽然已经观察到控制系统降低了振动的整体功率,但较高的分辨率谱表明,控制系统的另一个影响是在17.58 Hz的谐波分量陷波。因此,控制系统不仅降低了振动,而且使其更接近正弦。

值得注意的是,频率分辨率是由信号点数决定的,而不是由FFT点数决定的。增加FFT点的数量可以插值频率数据,为您提供更多频谱细节,但这并不能提高分辨率。

结论

在这个例子中,您学习了如何使用fftifft周期图pwelch,bandpower功能。你们理解了FFT的复杂性以及频谱的幅度和相位所包含的信息。在分析信号的周期性时,您看到了使用频域数据的优点。你们学习了如何计算一个噪声信号的总功率或特定频带的功率。您了解了如何提高频谱的频率分辨率使您能够观察间隔很近的频率成分,并且了解了频率分辨率和频谱平均之间的权衡。

进一步的阅读

有关频域分析的更多信息,请参阅信号处理工具箱。

参考文献:J.G. Proakis和d.g. Manolakis,“数字信号处理”。原理、算法与应用”,Prentice Hall, 1996。

附录

本例中使用了以下辅助函数。

另见

|||