Этот пример показывает, как получить непараметрические оценки спектральной плотности степени (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) аддитивном шуме. Синусоида имеет угловую частоту рад/образец. Используйте настройки по умолчанию генератора случайных чисел для воспроизводимых результатов.
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