exponenta event banner

Оценка спектральной плотности мощности с использованием БПФ

В этом примере показано, как получить оценки непараметрической спектральной плотности мощности (PSD), эквивалентные периодограмме с использованием fft. Примеры показывают, как правильно масштабировать выходные данные fft для входов четной длины, для нормированной частоты и герца и для односторонних и двусторонних оценок PSD.

Ввод четной длины с частотой выборки

Получение периодограммы для сигнала четной длины, дискретизированного на частоте 1 кГц, с использованием обоих fft и periodogram. Сравните результаты.

Создайте сигнал, состоящий из синусоидальной волны 100 Гц в N (0,1) аддитивном шуме. Частота дискретизации составляет 1 кГц. Длина сигнала составляет 1000 выборок. Используйте настройки по умолчанию генератора случайных чисел для воспроизводимых результатов.

rng default
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));

Получение периодограммы с помощью fft. Сигнал имеет действительное значение и четную длину. Поскольку сигнал является действительным, нужны только оценки мощности для положительных или отрицательных частот. Чтобы сохранить общую мощность, умножьте все частоты, возникающие в обоих множествах - положительную и отрицательную частоты - на коэффициент 2. Нулевая частота (DC) и частота Найквиста не встречаются дважды. Постройте график результата.

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

Figure contains an axes. The axes with title Periodogram Using FFT contains an object of type line.

Вычислите и постройте график периодограммы с помощью periodogram. Показать, что два результата идентичны.

periodogram(x,rectwin(length(x)),length(x),Fs)

Figure contains an axes. The axes with title Periodogram Power Spectral Density Estimate contains an object of type line.

mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x),Fs))
mxerr = 3.4694e-18

Вход с нормализованной частотой

Использовать fft получить периодограмму для ввода с использованием нормированной частоты. Создайте сигнал, состоящий из синусоидальной волны в N (0,1) аддитивном шуме. Синусоидальная волна имеет угловую частоту δ/4 рад/образец. Используйте настройки по умолчанию генератора случайных чисел для воспроизводимых результатов.

rng default
n = 0:999;
x = cos(pi/4*n) + randn(size(n));

Получение периодограммы с помощью fft. Сигнал имеет действительное значение и четную длину. Поскольку сигнал является действительным, нужны только оценки мощности для положительных или отрицательных частот. Чтобы сохранить общую мощность, умножьте все частоты, возникающие в обоих множествах - положительную и отрицательную частоты - на коэффициент 2. Нулевая частота (DC) и частота Найквиста не встречаются дважды. Постройте график результата.

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(2*pi*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:(2*pi)/N:pi;

plot(freq/pi,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Normalized Frequency (\times\pi rad/sample)') 
ylabel('Power/Frequency (dB/rad/sample)')

Figure contains an axes. The axes with title Periodogram Using FFT contains an object of type line.

Вычислите и постройте график периодограммы с помощью periodogram. Показать, что два результата идентичны.

periodogram(x,rectwin(length(x)),length(x))

Figure contains an axes. The axes with title Periodogram Power Spectral Density Estimate contains an object of type line.

mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x)))
mxerr = 1.4211e-14

Ввод комплексных значений с нормализованной частотой

Использовать fft получение периодограммы для комплексного входного сигнала с нормированной частотой. Сигнал является комплексным экспоненциальным с угловой частотой δ/4 рад/выборка в комплекснозначном N (0,1) шуме. Установите для генератора случайных чисел значения по умолчанию для воспроизводимых результатов.

rng default
n = 0:999;
x = exp(1j*pi/4*n) + [1 1j]*randn(2,length(n))/sqrt(2);

Использовать fft получить периодограмму. Так как входной сигнал является комплекснозначным, получайте периодограмму из [0,2δ) рад/выборки. Постройте график результата.

N = length(x);
xdft = fft(x);
psdx = (1/(2*pi*N)) * abs(xdft).^2;
freq = 0:(2*pi)/N:2*pi-(2*pi)/N;

plot(freq/pi,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Normalized Frequency (\times\pi rad/sample)') 
ylabel('Power/Frequency (dB/rad/sample)')

Figure contains an axes. The axes with title Periodogram Using FFT contains an object of type line.

Использовать periodogram для получения и построения графика периодограммы. Сравните оценки PSD.

periodogram(x,rectwin(length(x)),length(x),'twosided')

Figure contains an axes. The axes with title Periodogram Power Spectral Density Estimate contains an object of type line.

mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x),'twosided'))
mxerr = 2.8422e-14

См. также

Приложения

Функции