MATLAB高分辨率光谱分析

打开实时脚本

这个例子展示了如何在MATLAB®中使用高效滤波器组(有时被称为信道器)执行高分辨率频谱分析。为了进行比较，还展示了传统的平均修正周期图(Welch’s)方法。有关Simulink™中的类似示例，请参见在Simulink中的高分辨率光谱分析．

光谱分析的分辨率

在这种情况下，分辨率是指区分彼此附近的两个光谱成分的能力。分辨率取决于用于计算频谱的时域段的长度。当在时域段上使用窗口时，就像修改周期图的情况一样，所使用的窗口类型也会影响分辨率。

不同窗口的经典权衡是分辨率与旁瓣衰减之间的权衡。矩形窗口提供最高分辨率，但旁瓣衰减很差(~14 dB)。较差的旁瓣衰减会导致频谱分量被加窗操作所掩埋，因此是不可取的。汉窗以较低的频率分辨率为代价提供了良好的旁瓣衰减。Kaiser等可参数化窗口允许通过更改窗口参数来控制折衷。

而不是使用平均修正周期图(韦尔奇方法)，可以通过使用滤波器组方法来实现更高的分辨率估计，该方法模拟了模拟频谱分析仪的工作方式。其主要思想是使用滤波器组将信号分成不同的频率箱，并计算每个子带信号的平均功率。

基于滤波器组的频谱估计

对于本例，需要使用512个不同的带通滤波器来获得矩形窗口提供的相同分辨率。为了有效地实现512带通滤波器，使用了多相分析滤波器组(又名信道器)。它的工作原理是采用一个带宽为Fs/N的原型低通滤波器，其中N为所需的频率分辨率(本例中为512)，并以多相形式实现滤波器，就像实现FIR抽取器一样。不同于在十进制情况下将所有分支的结果相加，每个分支被用作n点FFT的输入。可以看出，FFT的每个输出对应一个低通滤波器的调制版本，从而实现了一个带通滤波器。滤波器组方法的主要缺点是由于多相滤波器增加了计算量，以及由于该滤波器的状态而对变化的信号的适应速度较慢。更多细节可以在fredric j. harris的《通信系统的多速率信号处理》一书中找到。Prentice Hall PTR, 2004。

在这个例子中，整个频谱估计使用了100个平均值。采样频率设置为1mhz。假设我们正在处理64个样本的帧，为了执行频谱估计，需要对这些帧进行缓冲。

NAvg = 100;Fs = 1e6;FrameSize = 64;NumFreqBins = 512;filterBankRBW = Fs/NumFreqBins;

简介实现了一个基于滤波器组的频谱估计器方法是相应设置的。在内部，它使用dsp。信道器实现多相滤波+ FFT(可用于频谱分析之外的其他应用，例如多载波通信)。

filterBankSA =频谱分析仪(.．.“SampleRate”Fs,.．.“RBWSource”，“属性”，.．.“RBW”filterBankRBW,.．.“AveragingMethod”，“指数”，.．.“ForgettingFactor”, 0.001,.．.“PlotAsTwoSidedSpectrum”假的,.．.“YLimits”50 [-150],.．.“YLabel”，“权力”，.．.“标题”，“滤波器组功率谱估计”，.．.“位置”，[50 375 800 450]);

测试信号

在本例中，测试信号是在64个样本帧中获取的。对于光谱分析而言，帧越大，分辨率越高。

测试信号由两个正弦波加上高斯白噪声组成。改变频率箱的数量，振幅，频率和噪声功率值是指导性的和鼓励的。

Sinegen = dsp。SineWave (“SampleRate”Fs,.．.“SamplesPerFrame”, FrameSize);

初始测试用例

首先，分别计算振幅为1和2的正弦波以及频率为200 kHz和250 kHz的滤波器组频谱估计。高斯白噪声的平均功率(方差)为1e-12。请注意，-114 dBm的单侧噪声底准确地显示在频谱估计中。

(sinegen) sinegen发布。振幅= [1 2];sinegen。频率= [200000 250000];noiseVar = 1e-12;% -114 dBm单侧noisfloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3);流(“噪声地板\ n”）;

噪声地板

流('滤波器组噪声底= %。2 f dBm \ n \ n”, noiseFloor);

滤波器组噪声底= -114.08 dBm

timesteps = 10 * ceil(NumFreqBins / FrameSize);为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);filterBankSA (x);结束发行版(filterBankSA)

利用谱估计器进行数值计算

dsp。SpectrumEstimator可用于计算滤波器组频谱估计。

为了给谱估计器提供更长的帧，在计算谱估计之前，缓冲区收集512个样本。虽然在本例中没有使用，但缓冲区允许重叠，这可以用于增加从给定数据集获得的平均值的数量。

filterBankEstimator = dsp。SpectrumEstimator (.．.“方法”，滤波器组的，.．.“AveragingMethod”，“指数”，.．.“ForgettingFactor”, 0.7,.．.“SampleRate”Fs,.．.“FrequencyRange”，“单向的”，.．.“PowerUnits”，dBm的）;buff = dsp.AsyncBuffer;release(sinegen) timesteps = 10 * ceil(NumFreqBins / FrameSize);为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);写(浅黄色,x);%缓冲数据如果迷。NumUnreadSamples >= NumFreqBins xbuff = read(buff,NumFreqBins);Pfbse = filterBankEstimator(xbuff);结束结束

比较使用不同方法的频谱估计

计算韦尔奇和滤波器组对振幅为1和2的正弦波和频率分别为200千赫和250千赫的频谱估计。高斯白噪声的平均功率(方差)为1e-12。

(sinegen) sinegen发布。振幅= [1 2];sinegen。频率= [200000 250000];filterBankSA。RBWSource =“汽车”；filterBankSA。健忘因子= 0.7;filterBankSA。位置= [50 375 400 450];光谱分析仪(.．.“方法”，“韦尔奇”，.．.“SampleRate”Fs,.．.“PlotAsTwoSidedSpectrum”假的,.．.“YLimits”50 [-150],.．.“YLabel”，“权力”，.．.“标题”，韦尔奇功率谱估计，.．.“位置”，[450 375 400 450]);noiseVar = 1e-12;timesteps = 500 * ceil(NumFreqBins / FrameSize);为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);filterBankSA (x);welchSA (x);结束

发行版(filterBankSA)

RBW = 488.28;hannNENBW = 1.5;welchNSamplesPerUpdate = Fs*hannNENBW/RBW;filterBankNSamplesPerUpdate = Fs/RBW;流(“样品/更新\ n”）;

样品/更新

流(韦尔奇样本/更新= %。3 f样品\ n”，.．.welchNSamplesPerUpdate);

Welch Samples/Update = 3072.008 Samples

流('过滤器组样本/更新= %。3 f样品\ n \ n”，.．.filterBankNSamplesPerUpdate);

Filter bank Samples/Update = 2048.005 Samples

welchnoisfloor = 10*log10((noiseVar/(welchNSamplesPerUpdate/2))/1e-3);filterbankisground = 10*log10((noiseVar/(filterBankNSamplesPerUpdate/2))/1e-3);流(“噪声地板\ n”）;

噪声地板

流(韦尔奇噪声地板= %。2 f dBm \ n”, welchNoiseFloor);

韦尔奇噪声底= -121.86 dBm

流('滤波器组噪声底= %。2 f dBm \ n \ n”, filterBankNoiseFloor);

滤波器组噪声底= -120.10 dBm

韦尔奇和基于滤波器组的频谱估计都在200千赫和250千赫检测到两种音调。基于滤波器组的频谱估计具有更好的音调隔离。对于相同的分辨率带宽(RBW)，平均修正周期图(Welch’s)方法需要3073个样本来计算频谱，而基于滤波器组的估计需要2048个样本。请注意，-120 dBm的单侧噪声底准确地显示在滤波器组光谱估计中。

比较使用不同窗口修改的周期图

考虑两种频谱分析仪，其中唯一的区别是使用的窗口:矩形或汉。

rectRBW = Fs/NumFreqBins;hannNENBW = 1.5;hannRBW = Fs*hannNENBW/NumFreqBins;矩形分析仪(.．.“方法”，“韦尔奇”，.．.“SampleRate”Fs,.．.“窗口”，“矩形”，.．.“RBWSource”，“属性”，.．.“RBW”rectRBW,.．.“PlotAsTwoSidedSpectrum”假的,.．.“YLimits”50 [-50],.．.“YLabel”，“权力”，.．.“标题”，基于矩形窗口的韦尔奇功率谱估计，.．.“位置”，[50 375 400 450]);光谱分析仪(.．.“方法”，“韦尔奇”，.．.“SampleRate”Fs,.．.“窗口”，“损害”，.．.“RBWSource”，“属性”，.．.“RBW”hannRBW,.．.“PlotAsTwoSidedSpectrum”假的,.．.“YLimits”50 [-150],.．.“YLabel”，“权力”，.．.“标题”，用汉窗估计韦尔奇功率谱，.．.“位置”，[450 375 400 450]);(sinegen) sinegen发布。振幅= [1 2];%也试试[0 2]sinegen。频率= [200000 250000];noiseVar = 1e-12;timesteps = 10 * ceil(NumFreqBins / FrameSize);为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);rectangularSA (x);hannSA (x);结束发行版(rectangularSA)

发行版(hannSA)

矩形窗口以低副瓣衰减为代价提供了窄的主瓣。相比之下，汉窗提供了更宽的主瓣，以换取更大的副瓣衰减。更宽的主瓣在250千赫时尤其明显。两个窗口都在正弦波所在的频率附近显示出较大的滚转。这可以掩盖低功率信号以上的噪声底。在滤波器组的情况下，这个问题几乎不存在。

将振幅更改为[0 2]而不是[1 2]有效地意味着有一个250 kHz的正弦波和噪声。这种情况很有趣，因为当200 kHz正弦波不干扰时，矩形窗口表现得特别好。原因是250khz是512个频率中的一个，平均划分1mhz。在这种情况下，由FFT中固有的频率采样引入的时域副本对用于功率谱计算的有时间限制的数据段进行了完美的周期扩展。一般来说，对于任意频率的正弦波，情况并非如此。这种对正弦波频率的依赖以及对信号干扰的敏感性是修正周期图方法的另一个缺点。

分辨率带宽(RBW)

一旦输入长度已知，每个分析仪的分辨率带宽就可以计算出来。RBW表示计算功率分量的带宽。也就是说，功率谱估计的每个元素代表了以估计的元素对应的频率为中心的宽度为RBW的带宽上的功率，单位为Watts、dBW或dBm。功率谱估计中每个元素的功率值是通过对跨频段的功率密度与RBW值进行积分得到的。RBW越低，分辨率越高，因为功率是在更细的网格(更小的带宽)上计算的。矩形窗口具有所有窗口中最高的分辨率。在凯撒窗的情况下，RBW取决于所使用的旁瓣衰减。

流(“RBW \ n”）

RBW

流(' welch - rectangle RBW = %。3 f赫兹\ n”, rectRBW);

韦尔奇矩形RBW = 1953.125 Hz

流('Welch-Hann RBW = %。3 f赫兹\ n”, hannRBW);

Welch-Hann RBW = 2929.688 Hz

流('滤波器组RBW = %。3 f赫兹\ n \ n”, filterBankRBW);

滤波器组RBW = 1953.125 Hz

在将振幅设置为[0 2]的情况下，即在250 kHz有一个正弦波的情况下，理解RBW和噪声下限之间的关系是很有趣的。预期的噪声底限为10*log10((noiseVar/(NumFreqBins/2))/1e-3)或约-114 dBm。矩形窗口对应的谱估计具有预期的噪声底，但使用汉窗的谱估计具有比预期高约2 dBm的噪声底。这是因为频谱估计是在512个频率点上计算的，但功率谱是在特定窗口的RBW上进行积分的。对于矩形窗口，RBW正好是1 MHz/512，因此频谱估计包含每个频率库的功率的独立估计。对于Hann窗口，RBW更大，因此频谱估计包含从一个频率库到下一个频率库的重叠功率。这种重叠功率增加了每个仓的功率值，并提高了噪声底限。可以分析计算如下:

hannnoiseflow = 10*log10((noiseVar/(NumFreqBins/2)*hannRBW/rectRBW)/1e-3);流(“噪声地板\ n”）;

噪声地板

流('汉噪声地板= %。2 f dBm \ n \ n”, hannNoiseFloor);

汉恩噪声底= -112.32 dBm

相互接近的正弦信号

为了说明解决问题，考虑下面的情况。正弦波频率改为200 kHz和205 kHz。滤波器组估计仍然准确。专注于基于窗口的估计器，由于Hann窗口的主瓣更宽，与矩形窗口估计相比，两个正弦信号更难区分。事实是，这两种估计都不是特别准确。请注意，205千赫基本上是我们能与200千赫区分的极限。对于更接近的频率，所有三个估计器都将无法分离两个频谱成分。分离较近的组件的唯一方法是有较大的帧尺寸，因此有较大的数量NumFrequencyBands在滤波器组估计器的情况下。

(sinegen) sinegen发布。振幅= [1 2];sinegen。频率= [200000 205000];filterBankSA。RBWSource =“属性”；filterBankSA。RBW= filterBankRBW; filterBankSA.Position = [850 375 400 450]; noiseVar = 1e-10; noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3);% -94 dBm单侧流(“噪声地板\ n”）;

噪声地板

流('噪音地板= %。2 f dBm \ n \ n”, noiseFloor);

噪声底= -94.08 dBm

timesteps = 500 * ceil(NumFreqBins / FrameSize);为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);filterBankSA (x);rectangularSA (x);hannSA (x);结束发行版(filterBankSA)

发行版(rectangularSA)

发行版(hannSA)

检测低功率正弦元件

接下来，重新运行前面的场景，但在170 kHz以非常小的振幅添加第三个正弦波。矩形窗估计和汉窗估计完全忽略了第三个正弦信号。滤波器组估计提供了更好的分辨率和更好的音调隔离，因此三个正弦波清晰可见。

(sinegen) sinegen发布。振幅= [1e-5 1 2];sinegen。频率= [170000 200000 205000];noiseVar = 1e-11;noisfloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3);% -104 dBm单侧流(“噪声地板\ n”）;

噪声地板

流('噪音地板= %。2 f dBm \ n \ n”, noiseFloor);

噪声底= -104.08 dBm

timesteps = 500 * ceil(NumFreqBins / FrameSize);为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);filterBankSA (x);rectangularSA (x);hannSA (x);结束发行版(filterBankSA)