В этом примере показано, как выполнить спектральный анализ высокого разрешения в MATLAB® с помощью эффективного набора фильтров, иногда называемого channelizer. В целях сравнения традиционная усредненная модифицированная периодограмма также показывают метод (валлийцев). Для подобного примера в Simulink™ смотрите Спектральный анализ Высокого разрешения в Simulink.
Разрешение в этом контексте относится к способности различать два спектральных компонента, которые лежат друг около друга. Разрешение зависит от длины сегмента временного интервала, использовался для расчета спектра. Когда работа с окнами используется на сегменте временного интервала, как имеет место с модифицированными периодограммами, тип окна, используемого также, влияет на разрешение.
Классическое согласование с различными окнами является одним из разрешения по сравнению с затуханием бокового лепестка. Прямоугольные окна обеспечивают самое высокое разрешение, но очень плохое затухание бокового лепестка (на ~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 который реализует многофазную фильтрацию плюс БПФ (и может использоваться для других приложений помимо анализа спектра, e.g. коммуникации мультинесущей).
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; % -114 dBm onesided noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3); fprintf('Noise Floor\n');
Noise Floor
fprintf('Filter bank noise floor = %.2f dBm\n\n',noiseFloor);Filter bank noise floor = -114.08 dBm
timesteps = 10 * ceil(NumFreqBins / FrameSize); for t = 1:timesteps x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1); filterBankSA(x); end release(filterBankSA)
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');Samples/Update
fprintf('Welch Samples/Update = %.3f Samples\n',... welchNSamplesPerUpdate);
Welch Samples/Update = 3072.008 Samples
fprintf('Filter bank Samples/Update = %.3f Samples\n\n',... filterBankNSamplesPerUpdate);
Filter bank Samples/Update = 2048.005 Samples
welchNoiseFloor = 10*log10((noiseVar/(welchNSamplesPerUpdate/2))/1e-3);
filterBankNoiseFloor = 10*log10((noiseVar/(filterBankNSamplesPerUpdate/2))/1e-3);
fprintf('Noise Floor\n');Noise Floor
fprintf('Welch noise floor = %.2f dBm\n',welchNoiseFloor);Welch noise floor = -121.86 dBm
fprintf('Filter bank noise floor = %.2f dBm\n\n',filterBankNoiseFloor);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)
Прямоугольное окно обеспечивает узкий основной лепесток за счет низкого затухания бокового лепестка. В отличие от этого окно Hann обеспечивает более широкий основной лепесток в обмен на намного большее затухание бокового лепестка. Более широкий основной лепесток особенно примечателен на уровне 250 кГц. Оба окна показывают большие спады вокруг частот, на которых лежат синусоиды. Это может замаскировать сигналы малой мощности интереса выше уровня шума. Та проблема фактически не существует в случае набора фильтров.
Изменение амплитуд к [0 2], а не [1 2] эффективно означает, что существует одна синусоида 250 кГц наряду с шумом. Этот случай интересен, потому что прямоугольное окно ведет себя особенно хорошо, когда синусоида на 200 кГц не вмешивается. Причина состоит в том, что 250 кГц являются одной из 512 частот, которые делят 1 МГц равномерно. В этом случае копии области времени, введенные частотой, производящей свойственный от БПФ, делают совершенное периодическое расширение сегмента ограниченных данных времени используемым для расчета спектра мощности. В общем случае для синусоид с произвольными частотами, дело обстоит не так. Эта зависимость от частоты синусоиды наряду с чувствительностью, чтобы сигнализировать об интерференции является другим недостатком модифицированного подхода периодограммы.
Полоса пропускания разрешения для каждого анализатора может быть вычислена, если входная длина известна. RBW указывает на полосу пропускания, по которой вычисляется компонент степени. То есть каждый элемент оценки спектра мощности представляет мощность в ваттах, dBW, или dBm через полосу пропускания ширины RBW, сосредоточенный вокруг частоты, соответствующей элементу оценки. Значение степени каждого элемента в оценке спектра мощности найдено путем интеграции плотности энергии по диапазону частот, заполненному значением RBW. Более низкий RBW указывает на более высокое разрешение, поскольку степень вычисляется по более прекрасной сетке (меньшая полоса пропускания). Прямоугольные окна имеют самое высокое разрешение всех окон. В случае окна Кайзера RBW зависит от используемого затухания бокового лепестка.
fprintf('RBW\n')RBW
fprintf('Welch-Rectangular RBW = %.3f Hz\n',rectRBW);Welch-Rectangular RBW = 1953.125 Hz
fprintf('Welch-Hann RBW = %.3f Hz\n',hannRBW);Welch-Hann RBW = 2929.688 Hz
fprintf('Filter bank RBW = %.3f Hz\n\n',filterBankRBW);Filter bank RBW = 1953.125 Hz
В случае установки амплитуд к [0 2], i.e. случай, в котором существует одна синусоида на уровне 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');Noise Floor
fprintf('Hann noise floor = %.2f dBm\n\n',hannNoiseFloor);Hann noise floor = -112.32 dBm
Чтобы проиллюстрировать проблему разрешения рассматривают следующий случай. Частоты синусоиды изменяются на 200 кГц и 205 кГц. Оценка набора фильтров остается точной. Фокусируясь на средствах оценки оконных, из-за более широкого основного лепестка в окне 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');
Noise Floor
fprintf('Noise floor = %.2f dBm\n\n',noiseFloor);Noise floor = -94.08 dBm
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)
Затем повторно выполните предыдущий сценарий, но добавьте третью синусоиду на уровне 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');
Noise Floor
fprintf('Noise floor = %.2f dBm\n\n',noiseFloor);Noise floor = -104.08 dBm
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)
