MATLAB中的高分辨率谱分析

打开实时脚本

此示例显示了如何在MATLAB®中使用高效滤波器组（有时称为通道化器）执行高分辨率光谱分析。为了便于比较，还显示了传统的平均修正周期图（Welch）方法。对于Simulink中的类似示例™, 看见万博1manbetxSimulink中的高分辨率光谱分析万博1manbetx．

光谱分析的分辨率

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

不同窗口的经典权衡是分辨率与旁瓣衰减。矩形窗口提供最高分辨率，但旁瓣衰减非常差（~14 dB）。旁瓣衰减差可能导致频谱成分被加窗操作掩埋，因此是不可取的。Hann窗提供了良好的旁瓣衰减，但牺牲了较低的频率分辨率。可参数化窗口（如Kaiser）允许通过更改窗口参数来控制权衡。

与使用平均修正周期图（Welch方法）不同，可以通过使用模拟频谱分析仪工作方式的滤波器组方法实现更高分辨率的估计。其主要思想是使用滤波器组将信号划分为不同的频率单元，并计算每个子带信号的平均功率。

基于滤波器组的频谱估计

对于这个例子，需要使用512个不同的带通滤波器来获得矩形窗口提供的相同分辨率。为了有效地实现512带通滤波器，使用了多相分析滤波器组（又称信道化器）。其工作原理是采用带宽为Fs/N的原型低通滤波器，其中N为所需的频率分辨率（本例中为512），并以多相形式实现滤波器，就像实现FIR抽取器一样。在抽取器的情况下，不是将所有分支的结果相加，而是将每个分支用作N点FFT的输入。可以看出，FFT的每个输出对应于低通滤波器的调制版本，因此实现了带通滤波器。滤波器组方法的主要缺点是由于多相滤波器而增加了计算量，以及由于滤波器的状态而对变化信号的适应较慢。更多细节可以在fredric j。哈里斯。Prentice Hall PTR，2004年。

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

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

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

filterBankSA = dsp。简介(．..“方法”，滤波器组的，．..“Numtapperband”, 24岁,．..“SampleRate”Fs,．..“RBWSource”，“属性”，．..“RBW”filterBankRBW,．..“SpectralAverages”NAvg,．..“绘制双侧面光谱”错误的．..“YLimits”50 [-150],．..“伊拉贝尔”，“权力”，．..“头衔”，“滤波器组功率谱估计”，．..“位置”，[50 375 800 450]);

测试信号

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

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

sinegen = dsp。SineWave (“SampleRate”Fs,．..“SamplesPerFrame”，框架尺寸）；

初始测试用例

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

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

噪声地板

流('滤波器组噪声下限= %。2 f dBm \ n \ n”，噪音地板）；

滤波器组噪声下限=-114.08 dBm

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

图频谱分析仪包含一个轴对象和其他类型的对象uiflowcontainer, uimenu, uitoolbar。标题为“滤波器组功率谱估计”的轴对象包含一个类型为line的对象。这个对象表示通道1。

使用频谱估计的数值计算

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

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

filterBankEstimator = dsp。SpectrumEstimator (．..“方法”，滤波器组的，．..“Numtapperband”, 24岁,．..“SampleRate”Fs,．..“SpectralAverages”NAvg,．..“频率范围”，“片面的”，．..“动力装置”，dBm的）;浅黄色= dsp.AsyncBuffer;释放(sinegen) timesteps = 10 * cell (NumFreqBins / FrameSize);为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);写(浅黄色,x);%缓冲数据如果buff.numreadsamples>=NumFreqBins xbuff=read（buff，NumFreqBins）；Pfbse=滤波器组刺激器（xbuff）；结束结束

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

计算振幅为1和2，频率为200 kHz和250 kHz的正弦波的welch和滤波器组谱估计。高斯白噪声的平均功率（方差）为1e-12。

(sinegen) sinegen发布。振幅= [1 2];sinegen。频率= [200000 250000];filterBankSA。RBWSource =“汽车”；filterBankSA.Position=[50 375 400 450]；welchSA=dsp.SpectrumAnalyzer(．..“方法”，“韦尔奇”，．..“SampleRate”Fs,．..“SpectralAverages”NAvg,．..“绘制双侧面光谱”错误的．..“YLimits”50 [-150],．..“伊拉贝尔”，“权力”，．..“头衔”，“韦尔奇功率谱估计”，．..“位置”，[450 375 400 450]);noiseVar = 1 e-12;timesteps = 500 * cell (NumFreqBins / FrameSize);为t=1:timesteps x=sum（sinegen（），2）+sqrt（noiseVar）*randn（FrameSize，1）；filterBankSA（x）；welchSA（x）；结束

Figure频谱分析仪包含一个轴对象和uiflowcontainer、uimenu、uitoolbar类型的其他对象。标题为Welch功率谱估计的axes对象包含一个line类型的对象。此对象表示通道1。

发行版(filterBankSA)

RBW = 488.28;hannNENBW = 1.5;welchNSamplesPerUpdate = f * hannNENBW / RBW;filterBankNSamplesPerUpdate = Fs / RBW;流('示例/更新\n'）;

样品/更新

流('Welch Samples/Update = %。3 f样品\ n”，．..welchNSamplesPerUpdate);

韦尔奇样本/更新=3072.008样本

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

过滤器组样本/更新= 2048.005样本

welchNoiseFloor=10*log10（（noiseVar/（welchNSamplesPerUpdate/2））/1e-3）；滤波器组噪声地板=10*log10（（噪声评估/（滤波器组样本更新/2））/1e-3）；fprintf(“噪声地板\ n”）;

噪声地板

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

Welch噪声底= -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；rectangularSA=dsp.SpectrumAnalyzer(．..“SampleRate”Fs,．..“窗口”，“矩形”，．..“RBWSource”，“属性”，．..“RBW”rectRBW,．..“SpectralAverages”NAvg,．..“绘制双侧面光谱”错误的．..“YLimits”,[-50 50],．..“伊拉贝尔”，“权力”，．..“头衔”，“使用矩形窗口的Welch功率谱估计”，．..“位置”,[50 375 400 450]); hannSA=dsp.SpectrumAnalyzer(．..“SampleRate”Fs,．..“窗口”，“损害”，．..“RBWSource”，“属性”，．..“RBW”，hannRBW，．..“SpectralAverages”NAvg,．..“绘制双侧面光谱”错误的．..“YLimits”50 [-150],．..“伊拉贝尔”，“权力”，．..“头衔”，“使用Hann窗的Welch功率谱估计”，．..“位置”，[450 375 400 450]);(sinegen) sinegen发布。振幅= [1 2];也试试[0 2]正弦频率=[200000 250000]；噪声评估=1e-12；时间步长=10*ceil（NumFreqBins/帧大小）；为t = 1:timesteps x = sum(sinegen()，2) + sqrt(noiseVar)*randn(FrameSize,1);rectangularSA (x);hannSA (x);结束发行版(rectangularSA)

图频谱分析仪包含一个轴对象和其他类型的对象uiflowcontainer, uimenu, uitoolbar。标题为Welch Power Spectrum Estimate using rectangle window的axis对象包含一个类型为line的对象。这个对象表示通道1。

释放（hannSA）

Figure Spectrum Analyzer包含axes对象和uiflowcontainer、uimenu、uitoolbar类型的其他对象。标题为Welch功率谱估计使用Hann窗口的axes对象包含line类型的对象。此对象表示通道1。

矩形窗口以低副瓣衰减为代价提供窄的主瓣。相比之下，Hann窗口提供更宽的主瓣以换取更大的副瓣衰减。在250 kHz时，更宽的主瓣尤其明显。两个窗口在正弦波l也就是说，这可以屏蔽噪声层以上的低功率信号。在滤波器组的情况下，这一问题实际上是不存在的。

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

分辨率带宽（RBW）

一旦输入长度已知，就可以计算每个分析仪的分辨率带宽。RBW表示计算功率分量的带宽。也就是说，功率谱估计的每个元素表示以功率为中心的带宽RBW上的功率（瓦特、dBW或dBm）估计的元素。功率谱估计中每个元素的功率值是通过积分RBW值跨越的频带上的功率密度来获得的。较低的RBW表示较高的分辨率，因为功率是在较细的网格（较小的带宽）上计算的.矩形窗口具有所有窗口中最高的分辨率。在凯撒窗口的情况下，RBW取决于使用的旁瓣衰减。

流(“RBW \ n”）

RBW

流('韦尔奇矩形RBW=%.3f Hz\n', rectRBW);

RBW = 1953.125 Hz

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

韦尔奇-汉恩RBW=2929.688赫兹

流('过滤银行RBW = %。3 f赫兹\ n \ n”，滤波器组（RBW）；

滤波器组RBW=1953.125 Hz

在将振幅设置为[0 2]的情况下，即在250 kHz处存在单一正弦波的情况下，了解RBW和噪声地板之间的关系是很有意思的。预期噪声地板为10*log10（（noiseVar/（NumFreqBins/2））/1e-3）或约-114 dBm。对应于矩形窗口的谱估计具有预期的噪声底，但使用Hann窗口的谱估计具有比预期高约2 dBm的噪声底。原因是谱估计在512个频率点计算，但功率谱在e特定窗口的RBW。对于矩形窗口，RBW正好为1 MHz/512，因此频谱估计包含每个频率单元的独立功率估计。对于Hann窗口，RBW更大，因此频谱估计包含从一个频率单元到下一个频率单元的重叠功率。此重叠功率增加设定每个箱子中的功率值，并提升噪声地板。该值可通过以下分析计算：

hannNoiseFloor=10*log10（（noiseVar/（NumFreqBins/2）*hannRBW/rectRBW）/1e-3）；fprintf(“噪声地板\ n”）;

噪声地板

流(“汉恩噪音地板= %。2 f dBm \ n \ n”, hannNoiseFloor);

汉恩噪音地板= -112.32 dBm

相互接近的正弦信号

为了说明解决问题，请考虑以下情况。正弦波频率分别为200khz和205khz。滤波器组的估计保持准确。聚焦于基于窗的估计，由于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”，噪音地板）；

噪声地板=-94.08 dBm

timesteps = 500 * cell (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处添加一个振幅非常小的第三个正弦波。矩形窗口估计和Hann窗口估计完全忽略了该第三个正弦波。滤波器组估计提供了更好的分辨率和更好的音调隔离，因此三个正弦波清晰可见。

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

噪声地板

流('噪音底限= %。2 f dBm \ n \ n”，噪音地板）；

噪音地板= -104.08 dBm

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