Этот пример показывает, как выполнять спектральный анализ с высоким разрешением с использованием эффективного набора фильтров, иногда называемого каналообразователем. Для целей сравнения также показан традиционный метод усредненной модифицированной периодограммы (Welch's).
Разрешение в этом контексте относится к способности различать две спектральные компоненты, которые лежат вблизи друг друга. Разрешение зависит от длины сегмента временной области, используемого для вычисления спектра. При использовании оконной обработки в сегменте временной области, как в случае с измененными периодограммами, тип используемого окна также влияет на разрешение.
Классический компромисс с различными окнами является одним из разрешения по сравнению с затуханием сайдлоба. Прямоугольные окна обеспечивают наивысшее разрешение, но очень низкое (~ 14 дБ) затухание боковых ребер. Слабое затухание боковых частей может привести к тому, что спектральные компоненты будут заглублены операцией оконной обработки и, таким образом, будут нежелательными. Окна Ганна обеспечивают хорошее затухание боковых зон за счет более низкого разрешения по частоте. Параметризуемые окна, такие как Kaiser, позволяют управлять компромиссом путем изменения параметра окна.
Вместо использования усредненных модифицированных периодограмм (метод Уэлча), более высокая оценка разрешения может быть достигнута с помощью подхода банка фильтров, который эмулирует работу аналоговых анализаторов спектра. Основная идея состоит в том, чтобы разделить сигнал на различные диапазоны частот, используя набор фильтров и вычисляя среднюю мощность каждого поддиапазонного сигнала.
Для этого примера необходимо использовать 512 различных полосовых фильтров для получения одинакового разрешения, обеспечиваемого прямоугольным окном. Чтобы эффективно реализовать 512 полосовых фильтров, набор фильтров многофазного анализа (a.k.a. каналообразователь). Это работает, принимая прототип фильтра нижних частот с шириной полосы частот Fs/N, где N желаемого разрешения частоты (512 в этом примере), и реализуя фильтр в многофазной форме, так же как реализован дециматор FIR. Вместо добавления результатов всех ветвей, как в случае прореживателя, каждая ветвь используется как вход в N-точечный БПФ. Можно показать, что каждый выход БПФ соответствует модулированной версии фильтра нижних частот, таким образом реализуя полосовой фильтр. Основным недостатком подхода к набору фильтров является увеличение вычислений за счет многофазного фильтра, а также более медленная адаптация к изменяющимся сигналам из-за состояний этого фильтра. Более подробную информацию можно найти в книге «Многоскоростная обработка сигналов для систем связи» fredric j. Харрис. Prentice Hall PTR, 2004.
В этом примере используются 100 средних значений спектральной оценки. Частота дискретизации устанавливается равной 1 МГц. Предполагается, что мы работаем с кадрами из 64 выборок, которые необходимо буферизовать для выполнения оценки спектра.
NAvg = 100; Fs = 1e6; FrameSize = 64; NumFreqBins = 512; filterBankRBW = Fs/NumFreqBins;
dsp.SpectrumAnalyzer реализует блок оценки спектра на основе банка фильтров, когда Method устанавливается соответствующим образом. Внутри, он использует dsp.Channelizer который реализует многофазную фильтрацию плюс FFT (и может использоваться для других приложений, помимо анализа спектра, например, связи с множеством несущих).
filterBankSA = dsp.SpectrumAnalyzer(... 'Method','Filter bank',... 'NumTapsPerBand',24,... 'SampleRate',Fs,... 'RBWSource','Property',... 'RBW',filterBankRBW,... 'SpectralAverages',NAvg,... 'PlotAsTwoSidedSpectrum',false,... 'YLimits',[-150 50],... 'YLabel','Power',... 'Title','Filter bank Power Spectrum Estimate',... 'Position',[50 375 800 450]);
В этом примере тестовый сигнал собирается в кадрах с 64 выборками. Для целей спектрального анализа чем больше кадр, тем лучше разрешение.
Тестовый сигнал состоит из двух синусоидальных волн плюс белый гауссов шум. Изменение количества частотных ячеек, амплитуды, частоты и значений мощности шума является поучительным и поощряемым.
sinegen = dsp.SineWave('SampleRate',Fs,'SamplesPerFrame',FrameSize);
Для начала вычисляют спектральную оценку банка фильтров для синусоидальных волн амплитуды 1 и 2 и частот 200 кГц и 250 кГц соответственно. Белый гауссов шум имеет среднюю мощность (дисперсию) 1e-12. Следует отметить, что в спектральной оценке точно показан единственный уровень шума -114 дБм.
release(sinegen) sinegen.Amplitude = [1 2]; sinegen.Frequency = [200000 250000]; noiseVar = 1e-12; noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3); % -114 dBm onesided fprintf('Noise Floor\n'); fprintf('Filter bank noise floor = %.2f dBm\n\n',noiseFloor); timesteps = 10 * ceil(NumFreqBins / FrameSize); for t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1); filterBankSA(x); end release(filterBankSA)
Noise Floor Filter bank noise floor = -114.08 dBm
dsp.SpectrumEstimator может использоваться для вычисления оценок спектра банка фильтров.
Для подачи в блок оценки спектра более длинного кадра буфер собирает 512 выборок перед вычислением спектральной оценки. Хотя в этом примере буфер не используется, он допускает перекрытие, которое может быть использовано для увеличения числа средних значений, полученных из данного набора данных.
filterBankEstimator = dsp.SpectrumEstimator(... 'Method','Filter bank',... 'NumTapsPerBand',24,... 'SampleRate',Fs,... 'SpectralAverages',NAvg,... 'FrequencyRange','onesided',... 'PowerUnits','dBm'); buff = dsp.AsyncBuffer; release(sinegen) timesteps = 10 * ceil(NumFreqBins / FrameSize); for t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1); write(buff,x); % Buffer data if buff.NumUnreadSamples >= NumFreqBins xbuff = read(buff,NumFreqBins); Pfbse = filterBankEstimator(xbuff); end end
Вычисляют спектральную оценку сварного шва и набора фильтров для синусоидальных волн амплитуды 1 и 2 и частот 200 кГц и 250 кГц соответственно. Белый гауссов шум имеет среднюю мощность (дисперсию) 1e-12.
release(sinegen) sinegen.Amplitude = [1 2]; sinegen.Frequency = [200000 250000]; filterBankSA.RBWSource = 'Auto'; filterBankSA.Position = [50 375 400 450]; welchSA = dsp.SpectrumAnalyzer(... 'Method','Welch',... 'SampleRate',Fs,... 'SpectralAverages',NAvg,... 'PlotAsTwoSidedSpectrum',false,... 'YLimits',[-150 50],... 'YLabel','Power',... 'Title','Welch Power Spectrum Estimate',... 'Position',[450 375 400 450]); noiseVar = 1e-12; timesteps = 500 * ceil(NumFreqBins / FrameSize); for t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1); filterBankSA(x); welchSA(x); end release(filterBankSA) RBW = 488.28; hannNENBW = 1.5; welchNSamplesPerUpdate = Fs*hannNENBW/RBW; filterBankNSamplesPerUpdate = Fs/RBW; fprintf('Samples/Update\n'); fprintf('Welch Samples/Update = %.3f Samples\n',welchNSamplesPerUpdate); fprintf('Filter bank Samples/Update = %.3f Samples\n\n',filterBankNSamplesPerUpdate); welchNoiseFloor = 10*log10((noiseVar/(welchNSamplesPerUpdate/2))/1e-3); filterBankNoiseFloor = 10*log10((noiseVar/(filterBankNSamplesPerUpdate/2))/1e-3); fprintf('Noise Floor\n'); fprintf('Welch noise floor = %.2f dBm\n',welchNoiseFloor); fprintf('Filter bank noise floor = %.2f dBm\n\n',filterBankNoiseFloor);
Samples/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
Оценка спектра на основе банка Welch и Filter обнаружила два тона с частотой 200 кГц и 250 кГц. Оценка спектра на основе банка фильтров имеет лучшую изоляцию тонов. Для того же самого способа разрешающей полосы пропускания (RBW) усредненная модифицированная периодограмма (Welch's) требует 3073 выборки для вычисления спектра по сравнению с 2048, требуемым для оценки банка фильтров. Следует отметить, что единичный уровень шума -120 дБм точно показан в спектральной оценке банка фильтров.
Рассмотрим два анализатора спектра, в которых единственным отличием является используемое окно: прямоугольное или ганновское.
rectRBW = Fs/NumFreqBins; hannNENBW = 1.5; hannRBW = Fs*hannNENBW/NumFreqBins; rectangularSA = dsp.SpectrumAnalyzer(... 'SampleRate',Fs,... 'Window','Rectangular',... 'RBWSource','Property',... 'RBW',rectRBW,... 'SpectralAverages',NAvg,... 'PlotAsTwoSidedSpectrum',false,... 'YLimits',[-50 50],... 'YLabel','Power',... 'Title','Welch Power Spectrum Estimate using Rectangular window',... 'Position',[50 375 400 450]); hannSA = dsp.SpectrumAnalyzer(... 'SampleRate',Fs,... 'Window','Hann',... 'RBWSource','Property',... 'RBW',hannRBW,... 'SpectralAverages',NAvg,... 'PlotAsTwoSidedSpectrum',false,... 'YLimits',[-150 50],... 'YLabel','Power',... 'Title','Welch Power Spectrum Estimate using Hann window',... 'Position',[450 375 400 450]); release(sinegen) sinegen.Amplitude = [1 2]; % Try [0 2] as well sinegen.Frequency = [200000 250000]; noiseVar = 1e-12; timesteps = 10 * ceil(NumFreqBins / FrameSize); for t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1); rectangularSA(x); hannSA(x); end release(rectangularSA) release(hannSA)
Прямоугольное окно обеспечивает узкий основной блок за счет низкого затухания боковых участков. Напротив, окно Ханна обеспечивает более широкий основной блок в обмен на гораздо большее затухание боковых блоков. Более широкий основной блок особенно заметен при 250 кГц. Оба окна демонстрируют большие откаты вокруг частот, на которых лежат синусоидальные волны. Это может маскировать сигналы низкой мощности, представляющие интерес, выше уровня шума. Эта проблема практически отсутствует в случае банка фильтров.
Изменение амплитуд на [0 2], а не [1 2] фактически означает наличие одиночной синусоидальной волны 250 кГц вместе с шумом. Этот случай интересен тем, что прямоугольное окно ведет себя особенно хорошо, когда синусоидальная волна 200 кГц не создает помех. Причина в том, что 250 кГц является одной из 512 частот, которые делятся 1 МГц равномерно. В этом случае реплики во временной области, вводимые частотной выборкой, присущей БПФ, делают идеальное периодическое расширение сегмента данных с ограниченным временем, используемого для вычисления спектра мощности. В общем случае для синусоидальных волн с произвольными частотами это не так. Эта зависимость от частоты синусоидальной волны наряду с восприимчивостью к помехам сигнала является еще одним недостатком модифицированного подхода к периодограмме.
Полоса пропускания разрешения для каждого анализатора может быть вычислена, как только длина входа известна. RBW указывает полосу пропускания, по которой вычисляется компонент мощности. То есть каждый элемент оценки спектра мощности представляет мощность в ваттах, дБВт или дБм по ширине полосы пропускания RBW, центрированной вокруг частоты, соответствующей элементу оценки. Значение мощности каждого элемента в оценке спектра мощности определяют путем интегрирования плотности мощности по полосе частот, охватываемой значением RBW. Более низкий RBW указывает на более высокое разрешение, поскольку мощность вычисляется по более тонкой сетке (меньшая полоса пропускания). Прямоугольные окна имеют наибольшее разрешение из всех окон. В случае окна Кайзера RBW зависит от используемого затухания боковых узлов.
fprintf('RBW\n') fprintf('Welch-Rectangular RBW = %.3f Hz\n',rectRBW); fprintf('Welch-Hann RBW = %.3f Hz\n',hannRBW); fprintf('Filter bank RBW = %.3f Hz\n\n',filterBankRBW);
RBW Welch-Rectangular RBW = 1953.125 Hz Welch-Hann RBW = 2929.688 Hz Filter bank RBW = 1953.125 Hz
В случае установки амплитуд на [0 2], т.е. в случае, когда имеется одна синусоидальная волна с частотой 250 кГц, интересно понять отношение между RBW и уровнем шума. Ожидаемый уровень шума равен 10 * log10 ((noureVar/( NumFreqBins/2 )/1e-3) или приблизительно -114 дБм. Спектральная оценка, соответствующая прямоугольному окну, имеет ожидаемый уровень шума, но спектральная оценка, использующая окно Ганна, имеет уровень шума, который примерно на 2 дБм выше ожидаемого. Причиной этого является то, что спектральная оценка вычисляется в 512 частотных точках, но спектр мощности интегрируется по RBW конкретного окна. Для прямоугольного окна RBW составляет ровно 1 MHz/512, так что спектральная оценка содержит независимые оценки мощности для каждого частотного бина. Для окна Ганна RBW больше, так что спектральная оценка содержит перекрывающуюся мощность от одного частотного бина к следующему. Эта перекрывающаяся мощность увеличивает значение мощности в каждом бункере и повышает уровень шума. Количество может быть вычислено аналитически следующим образом:
hannNoiseFloor = 10*log10((noiseVar/(NumFreqBins/2)*hannRBW/rectRBW)/1e-3); fprintf('Noise Floor\n'); fprintf('Hann noise floor = %.2f dBm\n\n',hannNoiseFloor);
Noise Floor Hann noise floor = -112.32 dBm
Для иллюстрации решения проблемы рассмотрим следующий случай. Частоты синусоиды меняют на 200 кГц и 205 кГц. Оценка банка фильтров остается точной. Фокусируясь на оконных оценках, из-за более широкого основного блока в окне Ганна, две синусоиды труднее различить по сравнению с прямоугольной оконной оценкой. Дело в том, что ни одна из двух оценок не является особенно точной. Обратите внимание, что 205 кГц в основном находится на пределе того, что мы можем отличить от 200 кГц. Для более близких частот все три оценщика не смогут разделить две спектральные компоненты. Единственный способ отделить более близкие компоненты - иметь большие размеры рамы и, следовательно, большее количество NumFrequencyBands в случае оценщика банка фильтров.
release(sinegen) sinegen.Amplitude = [1 2]; sinegen.Frequency = [200000 205000]; filterBankSA.RBWSource = 'Property'; filterBankSA.RBW = filterBankRBW; filterBankSA.Position = [850 375 400 450]; noiseVar = 1e-10; noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3); % -94 dBm onesided fprintf('Noise Floor\n'); fprintf('Noise floor = %.2f dBm\n\n',noiseFloor); timesteps = 500 * ceil(NumFreqBins / FrameSize); for t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1); filterBankSA(x); rectangularSA(x); hannSA(x); end release(filterBankSA) release(rectangularSA) release(hannSA)
Noise Floor Noise floor = -94.08 dBm
Затем повторно выполните предыдущий сценарий, но добавьте третью синусоиду при 170 кГц с очень малой амплитудой. Эта третья синусоида полностью пропускается оценкой прямоугольного окна и оценкой окна Ганна. Оценка банка фильтров обеспечивает как лучшее разрешение, так и лучшую изоляцию тонов, так что три синусоидальные волны хорошо видны.
release(sinegen) sinegen.Amplitude = [1e-5 1 2]; sinegen.Frequency = [170000 200000 205000]; noiseVar = 1e-11; noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3); % -104 dBm onesided fprintf('Noise Floor\n'); fprintf('Noise floor = %.2f dBm\n\n',noiseFloor); timesteps = 500 * ceil(NumFreqBins / FrameSize); for t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1); filterBankSA(x); rectangularSA(x); hannSA(x); end release(filterBankSA) release(rectangularSA) release(hannSA)
Noise Floor Noise floor = -104.08 dBm
Подобная настройка для спектральной оценки высокого разрешения, показанная выше, может быть смоделирована в Simulink с использованием блока Spectrum Analyzer. Модель «Фильтр-банк» иллюстрирует высокую разрешающую способность оценки спектра на основе банка фильтров и более низкий уровень шума по сравнению с методом Уэлча с использованием Simulink.
Рассмотрим следующий случай. Три синусоиды при 170 кГц, 200 кГц и 205kHz с амплитудами [1e-5 12]. Оценка прямоугольного окна полностью пропускает первую синусоиду. Оценка банка фильтров обеспечивает лучшее разрешение и лучшую изоляцию трех тонов.
open_system('SpectrumAnalyzerFilterBank'); sim('SpectrumAnalyzerFilterBank');
bdclose('SpectrumAnalyzerFilterBank');
Численные вычисления для спектральной оценки высокого разрешения, показанные выше, также могут быть смоделированы в Simulink с использованием блока спектральной оценки. Модель dspfilterbankspectrumessation иллюстрирует возможности высокого разрешения оценки спектра на основе банка фильтров и более низкий уровень шума по сравнению с методом Уэлча, используя Simulink. График массива используется для визуализации результатов. График массива обеспечивает удобный способ построения графика оценок спектра. Значения показаны в дБм, но вместо них можно легко использовать ватты или дБВт.
open_system('dspfilterbankspectrumestimation'); sim('dspfilterbankspectrumestimation');
bdclose('dspfilterbankspectrumestimation');