Оценка спектра мощности в MATLAB

Спектр мощности (PS) сигнала во временной области представляет собой распределение мощности, содержащейся в сигнале по частоте, на основе конечного набора данных. Представление сигнала в частотной области часто легче анализировать, чем представление во временной области. Многие приложения для обработки сигналов, такие как шумоподавление и идентификация системы, основаны на частотно-специфических модификациях сигналов. Целью спектральной оценки мощности является оценка спектра мощности сигнала из последовательности временных выборок. В зависимости от того, что известно о сигнале, методы оценки могут включать параметрические или непараметрические подходы и могут быть основаны на анализе во временной или частотной области. Например, общий параметрический метод включает подгонку наблюдений к авторегрессионной модели. Обычным непараметрическим методом является периодограмма. Спектр мощности оценивается с использованием методов преобразования Фурье, таких как метод Уэлча и метод банка фильтров. Для сигналов с относительно малой длиной подход набора фильтров дает спектральную оценку с более высоким разрешением, более точным уровнем шума и пиками более точными, чем метод Уэлча, с низкой спектральной утечкой или без нее. Эти преимущества достигаются за счет увеличения вычислений и более медленного отслеживания. Дополнительные сведения об этих методах см. в разделе Спектральный анализ. Можно также использовать другие методы, такие как метод максимальной энтропии.

В MATLAB ® можно выполнять спектральный анализ динамического сигнала в реальном времени с помощью dsp.SpectrumAnalyzer object™ системы. Спектральные данные можно просмотреть в анализаторе спектра и сохранить в переменной рабочей области с помощью isNewDataReady и getSpectrumData функции объекта. В качестве альтернативы можно использовать dsp.SpectrumEstimator Системный объект, за которым следует dsp.ArrayPlot объект для просмотра спектральных данных. Выходные данные dsp.SpectrumEstimator объект - спектральные данные. Эти данные могут быть получены для дальнейшей обработки.

Оцените спектр мощности с помощью dsp. SpectrumAnalyzer

Открыть сценарий

Для просмотра спектра мощности сигнала можно использовать dsp.SpectrumAnalyzer object™ системы. Вы можете изменить динамику входного сигнала и увидеть влияние этих изменений на спектр мощности сигнала в реальном времени.

Инициализация

Инициализируйте источник синусоидальной волны для генерации синусоидальной волны и анализатор спектра для отображения спектра мощности сигнала. Входная синусоидальная волна имеет две частоты: одну на 1000 Гц, а другую на 5000 Гц. Создать два dsp.SineWave объекты, один для генерации синусоидальной волны 1000 Гц, а другой для генерации синусоидальной волны 5000 Гц.

Fs = 44100;
Sineobject1 = dsp.SineWave('SamplesPerFrame',1024,'PhaseOffset',10,...
    'SampleRate',Fs,'Frequency',1000);
Sineobject2 = dsp.SineWave('SamplesPerFrame',1024,...
    'SampleRate',Fs,'Frequency',5000);
SA = dsp.SpectrumAnalyzer('SampleRate',Fs,'Method','Filter bank',...
    'SpectrumType','Power','PlotAsTwoSidedSpectrum',false,...
    'ChannelNames',{'Power spectrum of the input'},'YLimits',[-120 40],'ShowLegend',true);

Анализатор спектра использует подход набора фильтров для вычисления спектра мощности сигнала.

Оценка

Потоковая передача и оценка спектра мощности сигнала. Построить for-loop для выполнения 5000 итераций. В каждой итерации перетекают в 1024 выборки (один кадр) каждой синусоидальной волны и вычисляют спектр мощности каждого кадра. Для формирования входного сигнала добавьте две синусоидальные волны. Результирующий сигнал представляет собой синусоидальную волну с двумя частотами: одна на 1000 Гц, а другая на 5000 Гц. Добавьте шум Гаусса с нулевым средним значением и стандартным отклонением 0,001. Для получения спектральных данных для дальнейшей обработки используйте isNewDataReady и getSpectrumData функции объекта. Переменная data содержит спектральные данные, которые отображаются на анализаторе спектра вместе с дополнительными статистическими данными о спектре.

data = [];
for Iter = 1:7000
    Sinewave1 = Sineobject1();
    Sinewave2 = Sineobject2();
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    SA(NoisyInput);
     if SA.isNewDataReady
        data = [data;getSpectrumData(SA)];
     end
end
release(SA);

На выходе анализатора спектра можно увидеть два различных пика: один на 1000 Гц, а другой на 5000 Гц.

Разрешающая полоса пропускания (RBW) - это минимальная полоса частот, которая может быть разрешена анализатором спектра. По умолчанию RBWSource имущества dsp.SpectrumAnalyzer объект имеет значение Auto. В этом режиме RBW представляет собой отношение частотного диапазона к 1024. В двустороннем спектре это значение равно, в то $\frac{F_{s}}{1024}$ время как в одностороннем спектре оно равно. $\frac{\frac{F_{s}}{2}}{1024}$ Анализатор спектра в этом примере показывает односторонний спектр. Следовательно, RBW представляет собой (44100/2 )/1024 или 21.53Hz

Используя это значение, $RBW$ количество входных выборок, необходимых для вычисления одного спектрального обновления, $N_{samples}$ задается следующим уравнением: $N_{samples} = \frac{F_{s}}{RBW}$

В этом примере $N_{samples}$ используется 44100/21.53 или 2048 выборок.

$RBW$ вычисляется в режиме «Auto» и дает хорошее разрешение по частоте.

Чтобы различать две частоты в дисплее, расстояние между двумя частотами должно быть, по меньшей мере, RBW. В этом примере расстояние между двумя пиками составляет 4000 Гц, что больше, чем. $RBW$ Таким образом, вы можете видеть пики отчетливо. Измените частоту второй синусоидальной волны на 1015 Гц. Разница между этими двумя частотами меньше. $RBW$

release(Sineobject2);
Sineobject2.Frequency = 1015;
for Iter = 1:5000
    Sinewave1 = Sineobject1();
    Sinewave2 = Sineobject2();
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    SA(NoisyInput);
end
release(SA);

Вершины не различимы.

Для увеличения разрешения частоты уменьшите $RBW$ до 1 Гц.

SA.RBWSource = 'property';
SA.RBW = 1;
for Iter = 1:5000
    Sinewave1 = Sineobject1();
    Sinewave2 = Sineobject2();
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    SA(NoisyInput);
end
release(SA);

При увеличении теперь различимы два пика, которые находятся на расстоянии 15 Гц.

При увеличении разрешения частоты разрешение по времени уменьшается. Чтобы поддерживать хороший баланс между разрешением частоты и временным разрешением, измените RBWSource свойство для Auto.

Во время потоковой передачи можно изменить входные свойства или свойства анализатора спектра и немедленно увидеть эффект на выходе анализатора спектра. Например, измените частоту второй синусоидальной волны, когда индекс цикла кратен 1000.

release(Sineobject2);
SA.RBWSource = 'Auto';
for Iter = 1:5000
    Sinewave1 = Sineobject1();
    if (mod(Iter,1000) == 0)
        release(Sineobject2);
        Sineobject2.Frequency = Iter;
        Sinewave2 = Sineobject2();
    else
        Sinewave2 = Sineobject2();
    end
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    SA(NoisyInput);
end
release(SA);

При выполнении цикла потоковой передачи можно увидеть, что пик второй синусоидальной волны изменяется в соответствии со значением итерации. Аналогично, можно изменить любое из свойств анализатора спектра во время моделирования и увидеть соответствующее изменение в выходных данных.

Преобразование мощности между единицами измерения

Анализатор спектра содержит три блока для определения спектральной плотности мощности: Watts/Hz, dBm/Hz, и dBW/Hz. Соответствующими блоками питания являются Watts, dBm, и dBW. Для электротехнических приложений можно также просмотреть среднеквадратичное значение сигнала в Vrms или dBV. Тип спектра по умолчанию - Power in dBm.

Преобразование мощности в ваттах в дБВт и дБм

Мощность в dBW задается:

$_{PdBW} = 10log10 (мощность в ватт/1$ Вт)

Мощность в dBm задается:

$_{PdBm} = 10log10 (мощность в ватт/1$ милливатт)

Для синусоидального сигнала с амплитудой 1 В мощность одностороннего спектра в Watts задается:

$\begin{array}{l} _{PWatts} =^{} \\ _{} A2/2PWatts = \end{array}$ 1/2

В этом примере эта мощность равна 0,5 Вт. Соответствующая мощность в дБм задаётся:

$\begin{array}{l} _{PdBm} = 10log10 (мощность в ватт/1 \\ _{милливатт)} PdBm = {10log10}^{(} \end{array}$ 0,5/10 − 3)

Здесь мощность равна 26,9897 дБм. Чтобы подтвердить это значение с помощью поиска пиков, щелкните Инструменты > Измерения > Поиск пиков.

Для сигнала белого шума спектр является плоским для всех частот. Анализатор спектра в этом примере показывает односторонний спектр в диапазоне [0 Fs/2]. Для сигнала белого шума с дисперсией 1e-4 мощность на единицу полосы пропускания (_{Punitbandwidth}) равна 1e-4. Суммарная мощность белого шума в ваттах по всему диапазону частот задаётся:

$\begin{array}{l} _{Pwitenoise} =_{} Punitbandwidth * число \\ _{частотных ячеек,}^{Pwhitenoise} = (\frac{10 −}{4)} * \\ _{(Fs/2RBW),}^{} \frac{Pwhitenoise =}{(10} - \end{array}$ 4) * (2205021,53)

Количество частотных ячеек - это отношение общей полосы пропускания к RBW. Для одностороннего спектра общая полоса пропускания равна половине частоты дискретизации. RBW в этом примере составляет 21,53 Гц. При этих значениях суммарная мощность белого шума в ваттах равна 0,1024 Вт. В дБм мощность белого шума может быть вычислена с помощью 10 * log10 (0,1024/10 ^ -3), что равно 20,103 дБм.

Преобразование мощности в ваттах в dBFS

Если для спектральных единиц установлено значение dBFS и установить полную шкалу (FullScaleSourceКому Auto, мощность в dBFS вычисляется как:

$_{}_{PdBFS=20⋅log10} \sqrt{(_{}} Pwatts/Full_Scale$ )

где:

P_watts - мощность в ваттах
Для двойных и поплавковых сигналов Full_Scale является максимальным значением входного сигнала.
Для сигналов с фиксированной точкой или целочисленных сигналов Full_Scale является максимальным значением, которое может быть представлено.

Если указать полный масштаб вручную (набор FullScaleSource кому Property), питание в dBFS задается:

$_{}_{PFS=20⋅log10} \sqrt{(_{}} Pwatts/FS$ )

Где FS - полный коэффициент масштабирования, указанный в FullScale собственность.

Для синусоидального сигнала с амплитудой 1 В мощность одностороннего спектра в Watts задается:

$\begin{array}{l} _{PWatts} =^{} \\ _{} A2/2PWatts = \end{array}$ 1/2

В этом примере эта мощность равна 0,5 Вт, а максимальный входной сигнал для синусоидальной волны равен 1 В. Соответствующая мощность в дБФС задается следующим образом:

$_{}_{PFS=20⋅log10} \sqrt{(} 1/2/1$ )

Здесь мощность равна -3.0103. Чтобы подтвердить это значение в анализаторе спектра, выполните следующие команды:

Fs = 1000;  % Sampling frequency
sinef = dsp.SineWave('SampleRate',Fs,'SamplesPerFrame',100);
scope = dsp.SpectrumAnalyzer('SampleRate',Fs,...
   'SpectrumUnits','dBFS','PlotAsTwoSidedSpectrum',false)
%%
for ii = 1:100000
xsine = sinef();
scope(xsine)
end

Затем щелкните Инструменты > Измерения > Поиск пиков.

Преобразование мощности в дБм в среднеквадратичное значение в Vrms

Мощность в dBm задается:

$_{PdBm} = 10log10 (мощность в ватт/1$ милливатт)

Напряжение в СРК задается:

$_{Vrms} =^{_{}} \sqrt{^{10PdBm/2010}}$ − 3

В предыдущем примере _PdBm равен 26,9897 дБм. Значение _Vrms рассчитывается как

$_{Vrms} =^{} \sqrt{}$ 1026.9897/200.001

что равно 0,7071.

Для подтверждения этого значения:

Изменить тип на RMS.
Откройте поисковик пиков, щелкнув Инструменты > Измерения > Поисковик пиков.

Оцените спектр мощности с помощью dsp. SpectrumEstimator

Открыть сценарий

Альтернативно, можно вычислить спектр мощности сигнала, используя dsp.SpectrumEstimator Системный объект. Можно получить выходные данные оценщика спектра и сохранить данные для дальнейшей обработки. Просмотр других объектов в Estimation библиотека, тип help dsp в командной строке MATLAB ® и нажмитеEstimation.

Инициализация

Используйте тот же источник, что и в предыдущем разделе об использовании dsp.SpectrumAnalyzer для оценки спектра мощности. Входная синусоидальная волна имеет две частоты: одну на 1000 Гц, а другую на 5000 Гц. Инициализировать dsp.SpectrumEstimator вычисляют спектр мощности сигнала с использованием подхода набора фильтров. Просмотр спектра мощности сигнала с помощью dsp.ArrayPlot объект.

Fs = 44100;
Sineobject1 = dsp.SineWave('SamplesPerFrame',1024,'PhaseOffset',10,...
    'SampleRate',Fs,'Frequency',1000);
Sineobject2 = dsp.SineWave('SamplesPerFrame',1024,...
    'SampleRate',Fs,'Frequency',5000);

SpecEst = dsp.SpectrumEstimator('Method','Filter bank',...
    'PowerUnits','dBm','SampleRate',Fs,'FrequencyRange','onesided');
ArrPlot = dsp.ArrayPlot('PlotType','Line','ChannelNames',{'Power spectrum of the input'},...
    'YLimits',[-80 30],'XLabel','Number of samples per frame','YLabel',...
    'Power (dBm)','Title','One-sided power spectrum with respect to samples');

Оценка

Потоковая передача и оценка спектра мощности сигнала. Построить for-loop для выполнения 5000 итераций. В каждой итерации перетекают в 1024 выборки (один кадр) каждой синусоидальной волны и вычисляют спектр мощности каждого кадра. Добавьте гауссов шум со средним значением 0 и стандартным отклонением 0,001 к входному сигналу.

for Iter = 1:5000
    Sinewave1 = Sineobject1();
    Sinewave2 = Sineobject2();
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    PSoutput = SpecEst(NoisyInput);
    ArrPlot(PSoutput);
end

Используя подход набора фильтров, спектральная оценка имеет высокое разрешение, и пики являются точными без спектральной утечки.

Преобразовать ось X в частоту

По умолчанию график массива показывает спектральные данные мощности относительно количества выборок на кадр. Количество точек на оси X равно длине входного кадра. Анализатор спектра строит график спектральных данных мощности относительно частоты. Для одностороннего спектра частота изменяется в диапазоне [0 Fs/2]. Для двустороннего спектра частота изменяется в диапазоне [-Fs/2 Fs/2]. Чтобы преобразовать ось X графика массива из отсчетного в частотный, выполните следующие действия.

Щелкните значок «Свойства конфигурации».
Для одностороннего спектра - на вкладке «Главная» установите для параметра «Приращение образца» значение и смещение по оси X значение 0.
Для двустороннего спектра - на вкладке «Главная» установите для параметра «Приращение образца» значение, а для параметра «Смещение по оси X» значение.

В этом примере спектр является односторонним и, следовательно, приращение Sample и смещение X установлены на 44100/1024 и 0 соответственно. Чтобы задать частоту в кГц, установите для параметра Sample increment значение 44.1/1024.

ArrPlot.SampleIncrement = (Fs/1000)/1024;
ArrPlot.XLabel = 'Frequency (kHz)';
ArrPlot.Title = 'One-sided power spectrum with respect to frequency';

for Iter = 1:5000
    Sinewave1 = Sineobject1();
    Sinewave2 = Sineobject2();
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    PSoutput = SpecEst(NoisyInput);
    ArrPlot(PSoutput);
end

Обработка в реальном времени

Выходные данные dsp.SpectrumEstimator объект содержит спектральные данные и доступен для дальнейшей обработки. Данные могут обрабатываться в режиме реального времени или храниться в рабочей области.

Документация