主要内容

频域分析实用简介

这个例子展示了如何执行和解释基本的频域信号分析。该示例讨论了使用信号的频域表示和时域表示的优点,并使用模拟数据和真实数据说明了基本概念。这个例子回答了一些基本问题,例如:FFT的幅度和相位的含义是什么?信号是周期性的吗?我如何衡量权力?在这个波段中有一个或多个信号吗?

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

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

FFT的幅值和相位信息

信号的频域表示携带有关信号在每个频率上的幅度和相位的信息。这就是为什么FFT计算的输出是复杂的。一个复数,x美元,有一个实部,x_r美元,和虚部,x_i美元,以致于$x = x_r + jx_i$.的量级x美元计算为$ \√6 {(x_r ^ 2 + x_i ^ 2)} $的相位x美元计算为$ \反正切{(x_i / x_r)} $.你可以使用MATLAB函数腹肌而且分别得到任意复数的幅值和相位。

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

Fs = 44100;Y = audieread (“guitartune.wav”);

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

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

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

画出信号频谱的幅值和相位分量。震级以对数刻度(dB)绘制,很方便。属性来展开阶段打开函数,这样我们就能看到频率的连续函数。

magndey = abs(Y);% FFT的大小phaseY =展开(角度(Y));FFT的阶段helperFrequencyAnalysisPlot1 (F, magnitudeY phaseY NFFT)

你可以对频域向量Y应用傅里叶反变换来恢复时间信号。“对称”标志说明传输线你处理的是一个实数时间信号所以它会把逆变换中出现的虚分量归零因为计算中的数值不准确。请注意,原始时间信号y和恢复信号y1实际上是相同的(它们的差的范数是1e-14的数量级)。两者之间非常小的差异也是由于上面提到的数字不准确。播放并收听未转换的信号y1。

y1 = ifft(Y,NFFT,“对称”);规范(y-y1)
Ans = 3.7851e-14
hplayer = audioplayer(y1, Fs);玩(hplayer);

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

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

将滤波后的信号返回到时域传输线

ylp = ifft(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);

寻找信号周期

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

某办公楼温度循环特性分析

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

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

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

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

NFFT =长度(temp);% FFT点数F = (0: 1/NFFT: 1/2-1/NFFT)*F;频率向量TEMP = fft(TEMP,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²+ 0.5 v_i³$,在那里v_o美元输出电压和v_i美元为输入电压。以3.6 kHz的采样率捕获数据。输入v_i美元由一个单位振幅的60赫兹正弦波组成。由于非线性失真的性质,你应该期望放大器输出信号包含直流分量,60赫兹分量,以及120和180赫兹的第二和第三次谐波。

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

负载ampoutput1.matFs = 3600;NFFT = length(y);功率谱是在传递“功率”标志输入时计算的[P,F] =周期图(y,[],NFFT,Fs,“权力”);helperFrequencyAnalysisPlot2 (F, 10 * log10 (P),“以赫兹为单位的频率”...功率谱(dBW)[] [], [200 - -0.5])

功率谱图显示了DC、60和120 Hz的四个预期峰值中的三个。它还显示了更多的伪峰,这肯定是由信号中的噪声引起的。请注意,180赫兹的谐波完全被淹没在噪声中。

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

PdBW = 10*log10(P);power_at_DC_dBW = PdBW(F= 0)%瓦分贝[peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,“minpeakheight”, -11);peakFreqs_Hz = F(peakFreqIdx) peakPowers_dBW
power_at_DC_dBW = -7.8873 peakFreqs_Hz = 60 120 peakPowers_dBW = -0.3175 -10.2547

改进噪声信号的功率测量

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

加载一个由500e3个放大器输出点组成的更大的观测。保持用于执行fft的点数为3600,以便将地板(500e3/3600) = 138个fft求平均以获得功率谱。

负载ampoutput2.matSegmentLength = NFFT;功率谱是在传递“功率”标志输入时计算的[P,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,F,“权力”);helperFrequencyAnalysisPlot2 (F, 10 * log10 (P),“以赫兹为单位的频率”...功率谱(dBW)[] [], [200 - -0.5])

从情节上看,pwelch有效去除所有由噪声引起的杂散频率峰值。被噪声掩盖的180hz的光谱成分现在可见了。平均消除了频谱中的方差,这有效地产生更准确的功率测量。

测量一个频带的总平均功率和功率

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

PWR = sum(y.^2)/length(y)%瓦
PWR = 8.1697

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

pwr1 = sum(P)%瓦
Pwr1 = 8.1698

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

计算50hz ~ 70hz频段的功率。结果将包括60hz功率加上感兴趣波段的噪声功率:

pwr_band = bandpower(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,ones(SegmentLength,1),0,NFFT,F,psd的);pwr_band1 = bandpower(PSD,F,[50 70],psd的);pwr_band_dBW1 = 10*log10(pwr_band1)%瓦分贝
pwr_band_dBW1 = 0.0798

寻找光谱分量

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

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

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

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

开环和闭环加速度功率谱表明,控制系统有源时,加速度功率谱在4 ~ 11 dB之间减小。最大的衰减发生在约23.44千赫。降低11分贝意味着振动功率降低了12.6倍。总功率从0.1670瓦降低到0.059瓦,降低了2.83。

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

你在分析振动数据,你知道振动有循环行为。那么,为什么上面所示的频谱图不包含任何典型的循环行为的尖锐谱线呢?也许你错过了这些线,因为它们不能用64点线段长度获得的分辨率来解析?增加频率分辨率,看看是否有以前无法解析的谱线。中使用的数据段长度可以增加pwelch函数为512点。这产生了一个新的分辨率Fs/512 = 1.9531 Hz。在这种情况下,FFT平均值的数量减少到底部(10e3/512) = 19。显然,在使用时,平均值的数量和频率分辨率之间存在权衡pwelch.保持FFT点的数量为512。

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

请注意频率分辨率的提高是如何让您在开环频谱上观察到三个峰值,在闭环频谱上观察到两个峰值的。这些峰值以前是无法解决的。开环谱的峰间距约为11 Hz,小于长度为64的谱段的频率分辨率,但大于长度为512的谱段的分辨率。振动的循环行为现在是可见的。主振动频率为5.86 Hz,等距频率峰值表明它们是谐波相关的。虽然已经观察到控制系统降低了振动的整体功率,但更高分辨率的频谱表明,控制系统的另一个效果是在17.58 Hz处截断谐波成分。因此,控制系统不仅降低了振动,而且使其更接近正弦波。

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

结论

方法对信号进行频域分析fft传输线周期图pwelch,bandpower功能。你理解了FFT的复杂性质以及频谱的幅度和相位中包含的信息。在分析信号的周期性时,您看到了使用频域数据的优点。您学习了如何计算噪声信号的总功率或特定频带上的功率。您了解了提高频谱的频率分辨率如何使您能够观察到紧密间隔的频率成分,并了解了频率分辨率和频谱平均之间的权衡。

进一步的阅读

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

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

附录

本例中使用了以下helper函数。

另请参阅

|||