Main Content

FFT를 사용한 파워 스펙트럼 밀도 추정값

이 예제에서는 periodogram 함수와 fft 함수를 사용하여 동등한 비모수적 파워 스펙트럼 밀도(PSD) 추정값을 구하는 방법을 보여줍니다. 짝수 길이 입력값인 경우, 정규화 주파수와 헤르츠 단위 주파수인 경우, 단방향 및 양방향 PSD 추정값인 경우의 각기 다른 사례에서 fft의 출력값을 적절하게 스케일링하는 방법을 보여줍니다. 모든 사례에서 사각 윈도우를 사용합니다.

이 예제에서는 시간적으로 주파수 성분이 일정한 정상(stationary) 신호를 중점적으로 다룹니다. 시간 종속 주파수 성분을 갖는 비정상(nonstationary) 신호의 경우에는 단시간 푸리에 변환이 더 나은 분석 툴입니다. 자세한 내용은 Spectrogram Computation with Signal Processing Toolbox 항목을 참조하십시오.

샘플 레이트를 사용하는 짝수 길이 입력값

fftperiodogram을 모두 사용하여 1kHz로 샘플링된 짝수 길이 신호의 주기도를 구합니다. 결과를 비교합니다.

N(0,1) 가산성 잡음이 있는 100Hz 사인파로 구성된 신호를 생성합니다. 샘플링 주파수는 1kHz입니다. 신호 길이는 1000개 샘플입니다.

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,pow2db(psdx))
grid on
title("Periodogram Using FFT")
xlabel("Frequency (Hz)")
ylabel("Power/Frequency (dB/Hz)")

Figure contains an axes object. The axes object with title Periodogram Using FFT, xlabel Frequency (Hz), ylabel Power/Frequency (dB/Hz) contains an object of type line.

periodogram을 사용하여 주기도를 계산하고 플로팅합니다. 두 결과가 동일한지 표시합니다.

periodogram(x,rectwin(N),N,fs)

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate, xlabel Frequency (Hz), ylabel Power/frequency (dB/Hz) contains an object of type line.

mxerr = max(psdx'-periodogram(x,rectwin(N),N,fs))
mxerr = 3.4694e-18

정규화 주파수를 가지는 입력값

fft를 사용하여 정규화 주파수를 가지는 입력값에 대한 주기도를 생성합니다. N(0,1) 가산성 잡음이 있는 사인파로 구성된 신호를 생성합니다. 사인파는 π/4 rad/sample의 각주파수를 가집니다.

N = 1000;
n = 0:N-1;
x = cos(pi/4*n) + randn(size(n));

fft를 사용하여 주기도를 구합니다. 신호는 실수 값이며 짝수 길이를 가집니다. 신호가 실수 값이므로 양수 주파수나 음수 주파수에 대한 전력 추정값만 구하면 됩니다. 총 전력을 유지하려면 양의 주파수 세트와 음의 주파수 세트에 나타나는 모든 주파수에 인수 2를 곱하십시오. 영주파수(DC)와 나이퀴스트 주파수는 두 번 나타나지 않습니다. 결과를 플로팅합니다.

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,pow2db(psdx))
grid on
title("Periodogram Using FFT")
xlabel("Normalized Frequency (\times\pi rad/sample)")
ylabel("Power/Frequency (dB/(rad/sample))")

Figure contains an axes object. The axes object with title Periodogram Using FFT, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Power/Frequency (dB/(rad/sample)) contains an object of type line.

periodogram을 사용하여 주기도를 계산하고 플로팅합니다. 두 결과가 동일한지 표시합니다.

periodogram(x,rectwin(N),N)

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Power/frequency (dB/(rad/sample)) contains an object of type line.

mxerr = max(psdx'-periodogram(x,rectwin(N),N))
mxerr = 4.4409e-16

정규화 주파수를 가지는 복소수 값 입력 신호

fft를 사용하여 정규화 주파수를 가지는 복소수 값 입력 신호에 대한 주기도를 생성합니다. 이 신호는 복소수 값 N(0,1) 잡음을 가지며 π/4 rad/sample의 각주파수를 가지는 복소수 지수입니다.

N = 1000;
n = 0:N-1;
x = exp(1j*pi/4*n) + [1 1j]*randn(2,N)/sqrt(2);

fft를 사용하여 주기도를 구합니다. 입력 신호가 복소수 값이므로 [0,2π) rad/sample에서 주기도를 구합니다. 결과를 플로팅합니다.

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

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

Figure contains an axes object. The axes object with title Periodogram Using FFT, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Power/Frequency (dB/(rad/sample)) contains an object of type line.

periodogram을 사용하여 주기도를 구하고 플로팅합니다. PSD 추정값을 비교합니다.

periodogram(x,rectwin(N),N,"twosided")

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Power/frequency (dB/(rad/sample)) contains an object of type line.

mxerr = max(psdx'-periodogram(x,rectwin(N),N,"twosided"))
mxerr = 2.2204e-16

참고 항목

함수

외부 웹사이트