主要内容

高分辨率光谱分析

这个例子展示了如何使用一个高效的滤波器组来进行高分辨率的光谱分析,这个滤波器组有时被称为信道化器。为了便于比较,本文还展示了一种传统的平均修正周期图(Welch’s)方法。

频谱分析中的分辨率

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

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

与其使用平均修正周期图(Welch的方法),还可以使用一种模拟模拟频谱分析仪工作原理的滤波器组方法来实现更高分辨率的估计。其主要思想是用一个滤波器组将信号分成不同的频率箱,并计算每个子带信号的平均功率。

基于滤波器组的频谱估计

为此示例,512不同的带通滤波器需要用于获得由矩形窗口提供的相同分辨率。为了有效地实现512带通滤波器,使用多相分析滤波器组(A.k.a.Condoxizer)。这通过采用具有FS / N的带宽的原型低通滤波器,其中N个期望的频率分辨率(512在该示例中),并且实现了与FIR抽射器相同的多相形式的过滤器。而不是将所有分支的结果添加到Decimator案例中,而是每个分支用作N点FFT的输入。可以示出FFT的每个输出对应于低通滤波器的调制版本,从而实现带通滤波器。滤波器组方法的主要缺点是由于多相滤波器而增加的计算,以及对由于该滤波器的状态而改变信号的改变的适应性较慢。通过Fredric J的“关于通信系统”的“多速率信号处理”中可以找到更多细节。哈里斯。Prentice Hall Ptr,2004年。

在该示例中,整个频谱估计的100个平均值。采样频率设定为1 MHz。假设我们正在使用64个样本的帧,这需要缓冲以便执行频谱估计。

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

dsp。简介实现基于过滤器的频谱估计器方法相应的设置。在内部,它使用dsp.channelizer.这实现了多相滤波加FFT(并且除了频谱分析之外,可以用于其他应用,例如,多载波通信)。

filterbanksa = dsp.spectrumanalyzer(...“方法”滤波器组的...'numtapsperband', 24岁,...'采样率',fs,...“RBWSource”'财产'...'rbw',filterbankrbw,...'spectralaverages',Navg,...'plotastwosidedspectrum',错误的,...'ylimits'50 [-150],...'ylabel''力量'...'标题'“滤波器组功率谱估计”...'位置',[50 375 800 450]);

测试信号

在该示例中,在64样本帧中获取测试信号。对于光谱分析目的,框架越大,分辨率越好。

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

sinegen = dsp。SineWave ('采样率',fs,'samplesperframe',框架);

初始测试用例

要启动,分别计算幅度1和2的正弦波和200kHz和250kHz的正弦波频谱估计。白色高斯噪声具有1E-12的平均功率(方差)。请注意,在光谱估计中精确地显示-114 dBm的多个噪声底板。

释放(sinegen)sinegen.amplitude = [1 2];sinegen。频率= [200000 250000];noiseVar = 1 e-12;noiseFloor = 10 * log10 ((noiseVar / (NumFreqBins / 2)) / 1 e - 3);% -114 dBm单侧流('噪音地板\ n');流('过滤器银行噪声底板=%.2f dbm \ n \ n', noiseFloor);timesteps = 10 * cell (NumFreqBins / FrameSize);为了t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);Filterbanksa(x);结尾发行版(filterBankSA)
噪声层滤波器组噪声层= -114.08 dBm

使用频谱估计器的数值计算

dsp。SpectrumEstimator可以用来计算滤波器组的谱估计。

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

filterbankestimator = dsp.spectrumestimator(...“方法”滤波器组的...'numtapsperband', 24岁,...'采样率',fs,...'spectralaverages',Navg,...“FrequencyRange”'片面'...“PowerUnits”dBm的);浅黄色= dsp.AsyncBuffer;释放(sinegen) timesteps = 10 * cell (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 kHz和250 kHz的正弦波的welch和滤波器组谱估计值。白色高斯噪声具有1E-12的平均功率(方差)。

释放(sinegen)sinegen.amplitude = [1 2];sinegen。频率= [200000 250000];filterbanksa.rbwsource =.“汽车”;filterBankSA。Position = [50 375 400 450];welchSA = dsp。简介(...“方法”“韦尔奇”...'采样率',fs,...'spectralaverages',Navg,...'plotastwosidedspectrum',错误的,...'ylimits'50 [-150],...'ylabel''力量'...'标题'“韦尔奇功率谱估算”...'位置'[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);结尾释放(filterBankSA) RBW = 488.28;hannNENBW = 1.5;welchNSamplesPerUpdate = f * hannNENBW / RBW;filterBankNSamplesPerUpdate = Fs / RBW;流(“样品/更新\ n”);流('Welch Samples/Update = %。3 f样品\ n”,welchnsamplesperupdate);流('过滤器存储器样本/更新=%.3F样本\ n \ n',filterbanknsampleSperupdate);welchnoisefloor = 10 * log10((Noisevar /(WelchnsampleSperupdate / 2))/ 1E-3);filterbanknoisefloor = 10 * log10((Noisevar /(FilterbanknsampleSperupdate / 2))/ 1E-3);流('噪音地板\ n');流('韦尔奇噪音地板= %。2 f dBm \ n”,Welchnoisefloor);流('过滤器银行噪声底板=%.2f dbm \ n \ n', filterBankNoiseFloor);
sample /Update Welch Samples/Update = 3072.008 Samples Filter bank Samples/Update = 2048.005 Samples Noise Floor Welch Noise Floor = -121.86 dBm Filter bank Noise Floor = -120.10 dBm

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

使用不同的Windows比较修改的一段时间

考虑两种频谱分析仪,其中唯一不同的是使用的窗口:矩形或汉恩。

rectRBW = Fs / NumFreqBins;hannNENBW = 1.5;hannRBW = f * hannNENBW / NumFreqBins;rectangularSA = dsp。简介(...'采样率',fs,...'窗户'“矩形”...“RBWSource”'财产'...'rbw',recrbw,...'spectralaverages',Navg,...'plotastwosidedspectrum',错误的,...'ylimits'50 [-50],...'ylabel''力量'...'标题''Welch功率谱估计使用矩形窗口'...'位置',[50 375 400 450]);hannSA = dsp。简介(...'采样率',fs,...'窗户'“损害”...“RBWSource”'财产'...'rbw'hannRBW,...'spectralaverages',Navg,...'plotastwosidedspectrum',错误的,...'ylimits'50 [-150],...'ylabel''力量'...'标题'“利用Hann窗进行韦尔奇功率谱估计”...'位置'[450 375 400 450]);释放(sinegen)sinegen.amplitude = [1 2];也试试[0 2]sinegen。频率= [200000 250000];noiseVar = 1 e-12;timesteps = 10 * cell (NumFreqBins / FrameSize);为了t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);矩形(x);Hannsa(x);结尾释放(Rectangularsa)释放(Hannsa)

矩形窗口以降低旁瓣衰减为代价提供了狭窄的主瓣。相比之下,汉恩窗提供了更宽的主瓣,以换取更大的旁瓣衰减。更宽的主瓣在250千赫时特别明显。这两个窗口在正弦波所在的频率附近表现出很大的起伏。这可以掩盖低功率信号感兴趣的噪音地板以上。在滤波器组的情况下,这个问题几乎不存在。

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

分辨率带宽(RBW)

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

流('rbw \ n')流('Welch-矩形RBW =%.3F Hz \ n',recrbw);流('Welch-Hann RBW = %。3 f赫兹\ n”, hannRBW);流('过滤器银行RBW =%.3F Hz \ n \ n',Filterbankrbw);
RBW Welch-矩形RBW = 1953.125 Hz Welch-Hann RBW = 2929.688 Hz过滤器银行RBW = 1953.125 Hz

在将幅度设定为[02]的情况下,即在250 kHz处存在单个正弦波的情况,有趣的是了解RBW和噪声底板之间的关系。预期的噪声底板是10 * log10((Noisevar /(Numfreqbins / 2))/ 1E-3)或约-114 dBm。对应于矩形窗口的光谱估计具有预期的噪声底板,但使用HANN窗口的光谱估计具有比预期高约2 dBm的噪声底板。其原因是,光谱估计在512频率点计算,但功率谱集成在特定窗口的RBW上。对于矩形窗口,RBW精确为1 MHz / 512,因此光谱估计包含每个频率箱的功率的独立估计。对于HANN窗口,RBW更大,因此光谱估计包含从一个频率箱到下一个频率箱的重叠电力。该重叠功率增加了每个箱中的功率值并提升了噪声底板。该金额可以分析计算如下:

hannNoiseFloor = 10 * log10 ((noiseVar / (NumFreqBins / 2) * hannRBW / rectRBW) / 1 e - 3);流('噪音地板\ n');流(“汉恩噪音地板= %。2 f dBm \ n \ n”,hannnoisefloor);
噪声地板Hann噪音底线= -112.32 dBm

相互接近的正弦信号

为了说明解决问题,请考虑以下案例。正弦频率变为200 kHz和205 kHz。过滤器银行估计仍然是准确的。专注于基于窗口的估算器,因为HANN窗口中的大型阵容更广泛,与矩形窗口估计相比,两个正弦波更难区分。事实上,这两个估计都不是特别准确的。请注意,205 kHz基本上是我们可以区分200 kHz的限制。对于较近的频率,所有三个估算器都将无法分开两个光谱分量。分离较近组件的唯一方法是具有更大的帧大小,因此具有更大数量的NumFrequencyBands在滤波器组估计器的情况下。

释放(sinegen)sinegen.amplitude = [1 2];sinegen.frequency = [200000 205000];filterbanksa.rbwsource =.'财产';filterBankSA。RBW = filterBankRBW;filterBankSA。Position = [850 375 400 450];noiseVar = 1平台以及;noiseFloor = 10 * log10 ((noiseVar / (NumFreqBins / 2)) / 1 e - 3);%-94 dbm fonsided流('噪音地板\ n');流('噪音地板=%.2f dbm \ n', noiseFloor);timesteps = 500 * cell (NumFreqBins / FrameSize);为了t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);Filterbanksa(x);矩形(x);Hannsa(x);结尾发行版(filterBankSA)发布(rectangularSA)发布(hannSA)
噪声底噪声= -94.08 dBm

检测低功率正弦部件

接下来,重新运行前一个场景,但添加第三个正弦信号,频率为170 kHz,振幅很小。矩形窗估计和汉恩窗估计完全忽略了第三个正弦信号。滤波器组估计提供了更好的分辨率和更好的隔离音调,所以三个正弦波是清晰可见的。

(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');流('噪音地板=%.2f dbm \ n', noiseFloor);timesteps = 500 * cell (NumFreqBins / FrameSize);为了t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);Filterbanksa(x);矩形(x);Hannsa(x);结尾发行版(filterBankSA)发布(rectangularSA)发布(hannSA)
噪音地板噪音地板= -104.08 dBm

万博1manbetxSimulink版本的频谱分析仪

上面所示的高分辨率频谱估计的类似设置可以使用Simulink使用万博1manbetx频谱分析仪堵塞。这spectrumanalyzerfilterbank.模型说明了基于滤波器组的频谱估计的高分辨率能力和较低的噪声底限相比韦尔奇的方法,使用Simulink。万博1manbetx

考虑以下案例。具有170 kHz,200 kHz和205kHz的三个正弦曲线,具有振幅[1E-5 1 2]。第一个正弦曲线被矩形窗口估计完全错过。过滤器银行估计提供了更好的分辨率和更好地隔离三个音调。

示例模型

open_system (“SpectrumAnalyzerFilterBank”);SIM(“SpectrumAnalyzerFilterBank”);

关闭模型

bdclose(“SpectrumAnalyzerFilterBank”);

万博1manbetxSimulink版本的频谱估计器

上面所示的高分辨率频谱估计的数值计算也可以使用Simulink使用万博1manbetx谱估计堵塞。这dspfilterbankspectrumestimation.模型说明了基于滤波器组的频谱估计的高分辨率能力和较低的噪声底限相比韦尔奇的方法,使用Simulink。万博1manbetx数组图用于可视化结果。阵列图为谱估计的绘制提供了一种方便的方法。值以dBm表示,但瓦或dBW可以很容易地代替。

示例模型

open_system ('dspfilterbankspectrumestimation');SIM('dspfilterbankspectrumestimation');

关闭模型

bdclose('dspfilterbankspectrumestimation');