主要内容

频域分析的实用介绍

这个例子展示了如何执行和解释基本频域信号分析。该示例讨论了使用信号的频域表示和时域表示的优点,并使用模拟和真实数据说明了基本概念。这个例子回答了一些基本的问题,例如: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 = audioread(“guitartune.wav”);

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

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

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

画出信号频谱的幅值和相位分量。震级用对数刻度(dB)画出来很方便。方法将阶段解包装打开这样我们就能看到频率的连续函数。

magudey = 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输出中删除1khz以上的频率成分(使幅度等于零),并听这对音频文件的声音的影响。去除信号的高频成分称为低通滤波。

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

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

ylp = ifft(ylp,“对称”);

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

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)*Fs;频率矢量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^2 + 0.5 v_i^3$,在那里v_o美元输出电压和v_i美元为输入电压。采集数据的采样率为3.6 kHz。输入v_i美元由一个具有单位振幅的60赫兹正弦波组成。由于非线性失真的性质,你应该期望放大器输出信号包含直流分量,60hz分量,第二和第三次谐波在120和180 Hz。

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

负载ampoutput1.matFs = 3600;NFFT =长度(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,Fs,“权力”);helperFrequencyAnalysisPlot2 (F, 10 * log10 (P),“频率以赫兹为单位”...“功率谱(dBW)”[] [], [200 - -0.5])

正如图中所示,pwelch有效地去除所有由噪声引起的伪频率峰。在180赫兹的光谱成分被埋在噪音现在可见。平均从频谱中去除方差,这有效地产生更精确的功率测量。

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

测量时域信号的总平均功率是一项简单而常见的任务。对于放大器输出信号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,Fs,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,采样率为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分贝意味着振动功率降低了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赫兹”, {“开环”“闭环”}, 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,“数字信号处理”。原理、算法和应用”,普伦蒂斯霍尔,1996年。

附录

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

另请参阅

|||