В этом примере показано, как выполнить спектральный анализ высокого разрешения при помощи эффективного набора фильтров, иногда называемого channelizer. В целях сравнения традиционная усредненная модифицированная периодограмма также показывают метод (валлийцев).
Разрешение в этом контексте относится к способности различать два спектральных компонента, которые лежат друг около друга. Разрешение зависит от длины сегмента временного интервала, использовался для расчета спектра. Когда работа с окнами используется на сегменте временного интервала, как имеет место с модифицированными периодограммами, тип окна, используемого также, влияет на разрешение.
Классическое согласование с различными окнами является одним из разрешения по сравнению с затуханием бокового лепестка. Прямоугольные окна обеспечивают самое высокое разрешение, но очень плохое затухание бокового лепестка (на ~14 дБ). Плохое затухание бокового лепестка может привести к спектральным компонентам, прокладываемым под землей операцией работы с окнами, и таким образом является нежелательным. Окна Hann обеспечивают хорошее затухание бокового лепестка за счет более низкого разрешения частоты. Окна Parameterizable, такие как Кайзер позволяют управлять компромиссом путем изменения параметра окна.
Вместо того, чтобы использовать усредненные измененные периодограммы (метод валлийцев), более высокая оценка разрешения может быть достигнута при помощи подхода набора фильтров, который эмулирует, как аналоговый спектр анализаторы работает. Основная идея состоит в том, чтобы разделить сигнал на различные интервалы частоты с помощью набора фильтров и вычислив среднюю степень каждого сигнала поддиапазона.
В данном примере 512 различных полосовых фильтров должны использоваться, чтобы предоставить то же разрешение прямоугольным окном. Для того, чтобы реализовать эти 512 полосовых фильтров эффективно, многофазный аналитический набор фильтров (a.k.a. channelizer), используется. Это работает путем взятия прототипа фильтр lowpass с пропускной способностью Fs/N, где N желаемое разрешение частоты (512 в этом примере), и реализация фильтра в многофазной форме во многом как КИХ decimator реализован. Вместо того, чтобы добавить результаты всех ветвей как в decimator случае, каждая ветвь используется в качестве входа к БПФ N значений. Можно показать, что каждый выход БПФ соответствует модулируемая версия фильтра lowpass, таким образом реализовывая полосовой фильтр. Основной недостаток подхода набора фильтров является увеличенным расчетом из-за многофазного фильтра, а также более медленной адаптации к изменению сигналов из-за состояний того фильтра. Больше деталей может быть найдено в книге 'Многоскоростной Обработкой сигналов для Систем связи' fredric j. harris. PTR 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 dBm точно показывают в спектральной оценке.
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
Обе валлийских и основанных на наборе фильтров оценки спектра обнаружили два тона на уровне 200 кГц и 250 кГц. Основанная на наборе фильтров оценка спектра имеет лучшую изоляцию тонов. Для той же Пропускной способности Разрешения (RBW) усредненная измененная периодограмма метод (валлийцев) требует, чтобы 3 073 выборки вычислили спектр по сравнению с 2 048 требуемыми основанной на наборе фильтров оценкой. Обратите внимание на то, что односторонний уровень шума-120 dBm точно показывают в наборе фильтров спектральную оценку.
Считайте два спектра анализаторами, в которых единственной разницей является используемое окно: прямоугольный или 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)
Прямоугольное окно обеспечивает узкий mainlobe за счет низкого затухания бокового лепестка. В отличие от этого окно Hann обеспечивает более широкий mainlobe в обмен на намного большее затухание бокового лепестка. Более широкий mainlobe особенно примечателен на уровне 250 кГц. Оба окна показывают большие спады вокруг частот, на которых лежат синусоиды. Это может замаскировать сигналы малой мощности интереса выше уровня шума. Та проблема фактически не существует в случае набора фильтров.
Изменение амплитуд к [0 2], а не [1 2] эффективно означает, что существует одна синусоида 250 кГц наряду с шумом. Этот случай интересен, потому что прямоугольное окно ведет себя особенно хорошо, когда синусоида на 200 кГц не вмешивается. Причина состоит в том, что 250 кГц являются одной из 512 частот, которые делят 1 МГц равномерно. В этом случае копии области времени, введенные частотой, производящей свойственный от БПФ, делают совершенное периодическое расширение сегмента ограниченных данных времени используемым в расчете спектра мощности. В общем случае для синусоид с произвольными частотами, дело обстоит не так. Эта зависимость от частоты синусоиды наряду с чувствительностью, чтобы сигнализировать об интерференции является другим недостатком модифицированного подхода периодограммы.
Пропускная способность разрешения для каждого анализатора может быть вычислена, если входная длина известна. RBW указывает на пропускную способность, по которой вычисляется компонент степени. То есть каждый элемент оценки спектра мощности представляет мощность в ваттах, dBW, или dBm через пропускную способность ширины 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 dBm. Спектральная оценка, соответствующая прямоугольному окну, имеет ожидаемый уровень шума, но спектральная оценка с помощью окна Hann имеет уровень шума, который является приблизительно 2 dBm выше, чем ожидалось. Причина этого состоит в том, что спектральная оценка, вычисляют в 512 точках частоты, но спектр мощности интегрирован по RBW конкретного окна. Для прямоугольного окна RBW - точно 1 МГц / 512 так, чтобы спектральная оценка содержала независимые оценки степени для каждого интервала частоты. Для окна Hann 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 кГц. Оценка набора фильтров остается точной. Фокусируясь на средствах оценки оконных, из-за более широкого mainlobe в окне Hann, эти две синусоиды более трудно отличить, когда по сравнению с прямоугольным окном оценивают. Факт - то, что ни одна из двух оценок не особенно точна. Обратите внимание на то, что 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 кГц с очень маленькой амплитудой. Эта третья синусоида пропущена полностью прямоугольной оценкой окна и оценкой окна Hann. Оценка набора фильтров обеспечивает и лучшее разрешение и лучшую изоляцию тонов, так, чтобы эти три синусоиды ясно отобразились.
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 кГц и 205 кГц с амплитудами [1e-5 1 2]. Первая синусоида полностью пропущена прямоугольной оценкой окна. Оценка набора фильтров обеспечивает лучшее разрешение и лучшую изоляцию трех тонов.
open_system('SpectrumAnalyzerFilterBank'); sim('SpectrumAnalyzerFilterBank');
bdclose('SpectrumAnalyzerFilterBank');
Численные расчеты для высокого разрешения спектральная оценка, показанная выше, могут также быть смоделированы в Simulink с помощью блока Spectrum Estimator. dspfilterbankspectrumestimation модель иллюстрирует высокие разрешающие способности основанной на наборе фильтров оценки спектра и более низкого уровня шума по сравнению с методом валлийцев, с помощью Simulink. График массивов используется, чтобы визуализировать результаты. График массивов обеспечивает удобный способ построить оценки спектра. Значения показывают в dBm, но Уоттс или dBW могли легко использоваться вместо этого.
open_system('dspfilterbankspectrumestimation'); sim('dspfilterbankspectrumestimation');
bdclose('dspfilterbankspectrumestimation');