В этом примере показано, как получить оценки непараметрической спектральной плотности мощности (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)')

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

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)')

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

mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x)))
mxerr = 1.4211e-14
Использовать fft получение периодограммы для комплексного входного сигнала с нормированной частотой. Сигнал является комплексным экспоненциальным с угловой частотой рад/выборка в комплекснозначном N (0,1) шуме. Установите для генератора случайных чисел значения по умолчанию для воспроизводимых результатов.
rng default
n = 0:999;
x = exp(1j*pi/4*n) + [1 1j]*randn(2,length(n))/sqrt(2);Использовать fft получить периодограмму. Так как входной сигнал является комплекснозначным, получайте периодограмму из рад/выборки. Постройте график результата.
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)')

Использовать periodogram для получения и построения графика периодограммы. Сравните оценки PSD.
periodogram(x,rectwin(length(x)),length(x),'twosided')
mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x),'twosided'))mxerr = 2.8422e-14
fft | periodogram | pspectrum