Этот пример показывает, как выполнить спектральный анализ с высоким разрешением с использованием эффективной группы фильтров, иногда называемой канализатором. Для целей сравнения также показан традиционный усредненный модифицированный метод периодограммы (Welch's).
Разрешение в этом контексте относится к способности различать два спектральных компонента, которые находятся поблизости друг от друга. Разрешение зависит от длины сегмента временной области, используемого для вычисления спектра. Когда оконная обработка используется на сегменте временной области, как в случае с измененными периодограммами, тип используемого окна также влияет на разрешение.
Классический компромисс с различными окнами является одним из разрешения по сравнению с ослаблением бокового колеса. Прямоугольные окна обеспечивают самое высокое разрешение, но очень плохое (~ 14 дБ) ослабление бокового колеса. Плохое ослабление бокового elobe может привести к тому, что спектральные компоненты будут зарыты операцией оконной обработки и, таким образом, нежелательны. Окна Hann обеспечивают хорошее ослабление бокового колеса за счет более низкого разрешения частоты. Параметризуемые окна, такие как Kaiser, позволяют управлять компромиссом путем изменения параметра окна.
Вместо использования усредненных модифицированных периодограмм (метод Уэлча), более высокая оценка разрешения может быть достигнута с помощью подхода банка фильтров, который эмулирует работу аналоговых анализаторов спектра. Основная идея состоит в том, чтобы разделить сигнал на различные интервалы частот с помощью набора фильтров и вычисления средней степени каждого поддиапазона сигнала.
В данном примере необходимо использовать 512 различных полосно-пропускающих фильтров, чтобы получить то же разрешение, которое дает прямоугольное окно. В порядок эффективной реализации 512 полосно-пропускающих фильтров, набор фильтров анализа полифазы (a.k.a. channelizer). Это работает путем взятия прототипа lowpass фильтра с шириной полосы Fs/N, где N желаемое разрешение частоты (512 в этом примере), и реализации фильтра в полифазной форме, во многом подобной конечной импульсной характеристике. Вместо добавления результатов всех ветвей, как в случае дециматора, каждая ветвь используется как вход к N значений. Можно показать, что каждый выход БПФ соответствует модулированной версии lowpass, таким образом реализуя полосно-пропускающий фильтр. Основным недостатком подхода группы фильтров является увеличение расчетов из-за полифазного фильтра, а также более медленная адаптация к изменяющимся сигналам из-за состояний этого фильтра. Более подробную информацию можно найти в книге 'Multirate Signal Processing for Коммуникационные Системы' по fredic j. Харрис. Prentice Hall ПТР, 2004.
В этом примере используются 100 средних значений оценки спектра. Частота дискретизации устанавливается на 1 МГц. Принято, что мы работаем с системами координат из 64 выборок, которые необходимо будет буферизировать порядок для выполнения оценки спектра.
NAvg = 100; Fs = 1e6; FrameSize = 64; NumFreqBins = 512; filterBankRBW = Fs/NumFreqBins;
dsp.SpectrumAnalyzer реализует оценку спектра на основе банка фильтров при Method устанавливается соответствующим образом. Внутренне он использует dsp.Channelizer который реализует фильтрацию полифазы плюс БПФ (и может использоваться в других приложениях, помимо спектра анализа, например, связи с несколькими несущими).
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 Bank обнаружила два тональных сигнала на частотах 200 кГц и 250 кГц. Оценка спектра на основе банка фильтров имеет лучшую изоляцию тонов. Для той же полосы пропускания разрешения (RBW) усредненный модифицированный метод периодограммы (Welch's) требует 3073 выборки для вычисления спектра по сравнению с 2048, требуемым оценкой на основе банка фильтров. Обратите внимание, что одноуровневый шумовой пол -120 дБм точно показан в спектральной оценке группы фильтров.
Рассмотрим два спектра анализатора, в которых единственным различием является используемое окно: прямоугольный или Hann.
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 ((noiseVar/( 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. Модель SpectrumAnalyzerFilterBank иллюстрирует высокие возможности разрешения при оценке спектра на основе фильтра и более низкий шумовой пол по сравнению с методом Уэлча, используя Simulink.
Рассмотрим следующий случай. Три синусоиды при 170 кГц, 200 кГц и 205kHz с амплитудами [1e-5 1 2]. Первая синусоида полностью пропущена прямоугольной оценкой окна. Оценка банка фильтров обеспечивает лучшее разрешение и лучшую изоляцию трех тонов.
open_system('SpectrumAnalyzerFilterBank'); sim('SpectrumAnalyzerFilterBank');
bdclose('SpectrumAnalyzerFilterBank');
Численные расчеты для спектральной оценки высокого разрешения, показанные выше, могут также быть смоделированы в Simulink с использованием блока Spectrum Estimator. Модель dspfilterbankspectrumestimation иллюстрирует возможности высокого разрешения оценки спектра на основе банка фильтров и более низкого уровня шума по сравнению с методом Уэлча, используя Simulink. График массива используется для визуализации результатов. График массива обеспечивает удобный способ построения графиков оценок спектра. Значения показаны в дБм, но вместо этого можно легко использовать Watts или dBW.
open_system('dspfilterbankspectrumestimation'); sim('dspfilterbankspectrumestimation');
bdclose('dspfilterbankspectrumestimation');