В этом примере показано, как получить непараметрические оценки спектральной плотности мощности (PSD), эквивалентные периодограмме с помощью fft
. Примеры показывают вам, как правильно масштабировать выход fft
для входных параметров ровной длины, для нормированной частоты и герц, и для одного - и двухсторонние оценки PSD.
Получите периодограмму для сигнала ровной длины, произведенного на уровне 1 кГц с помощью обоих fft
и periodogram
Сравнение результатов.
Создайте сигнал, состоящий из синусоиды на 100 Гц в N (0,1) аддитивный шум. Частота дискретизации составляет 1 кГц. Длина сигнала является 1 000 выборок. Используйте настройки по умолчанию генератора случайных чисел для восстанавливаемых результатов.
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) аддитивный шум. Синусоида имеет угловую частоту рад/отсчет. Используйте настройки по умолчанию генератора случайных чисел для восстанавливаемых результатов.
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