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

Этот пример показывает, как получить непараметрические оценки спектральной плотности степени (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

См. также

Приложения

Функции