periodogram
주기도 파워 스펙트럼 밀도 추정값
구문
설명
[는 위에 열거된 구문에 주파수 범위와 스펙트럼 유형도 지정합니다. 이러한 입력 인수 중 하나를 지정하거나 둘 다 지정할 수 있습니다.pxx,f] = periodogram(___,freqRange,spectrumType)
periodogram(___)(출력 인수 없음)은 현재 Figure 창에 주기도 PSD 추정값 또는 파워 스펙트럼을 플로팅합니다.
예제
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 주기도를 구합니다.
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 사인파를 생성합니다. 신호의 길이는 320개 샘플입니다. 디폴트 사각 윈도우와 DFT 길이를 사용하여 주기도를 구합니다. DFT 길이는 신호 길이보다 큰 첫 번째 2의 거듭제곱값이거나 512개 점입니다. 신호가 실수 값을 가지고 짝수 길이이기 때문에 주기도가 단측이고 512/2+1개 점이 있습니다.
n = 0:319; x = cos(pi/4*n) + randn(size(n)); [pxx,w] = periodogram(x); plot(w/pi,pow2db(pxx)) xlabel("Normalized Frequency (\times \pi rad/sample)") title("Periodogram Power Spectral Density Estimate")
출력값을 지정하지 않고 periodogram을 사용하여 다시 플롯을 그립니다.
periodogram(x)
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 수정된 주기도를 구합니다.
가산성 백색 잡음과 함께 각주파수 radian/sample을 갖는 사인파를 생성합니다. 신호의 길이는 320개 샘플입니다. 해밍 윈도우와 디폴트 DFT 길이를 사용하여 수정된 주기도를 구합니다. DFT 길이는 신호 길이보다 큰 첫 번째 2의 거듭제곱값이거나 512개 점입니다. 신호가 실수 값을 가지고 짝수 길이이기 때문에 주기도가 단측이고 512/2+1개 점이 있습니다.
n = 0:319; x = cos(pi/4*n) + randn(size(n)); periodogram(x,hamming(length(x)))
가산성 백색 잡음과 함께 각주파수 radian/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 주기도를 구합니다. 신호 길이와 동일한 DFT 길이를 사용합니다.
가산성 백색 잡음과 함께 각주파수 radian/sample을 갖는 사인파를 생성합니다. 신호의 길이는 320개 샘플입니다. 디폴트 사각 윈도우 및 신호 길이와 같은 DFT 길이를 사용하여 주기도를 구합니다. 신호가 실수 값을 가지기 때문에 기본적으로 길이가 320/2+1인 단측 주기도가 반환됩니다.
rng("default")
n = 0:319;
x = cos(pi/4*n) + randn(size(n));
nfft = length(x);
periodogram(x,[],nfft)
1700년부터 1987년까지 매년 샘플링된 볼프(상대 흑점)수 데이터에 대한 주기도를 구합니다.
상대 흑점수 데이터를 불러옵니다. 디폴트 사각 윈도우와 DFT 점 개수(이 예제의 경우 512)를 사용하여 주기도를 구합니다. 이 데이터의 샘플 레이트는 1샘플/년입니다.
load sunspot.dat
relNums=sunspot(:,2);
[pxx,f] = periodogram(relNums,[],[],1);주기도를 플로팅합니다. 주기도의 약 0.1주기/년 지점에 피크가 있으며, 이는 대략 10년의 기간을 나타냅니다.
plot(f,pow2db(pxx)) xlabel("Cycles/Year") ylabel("dB / (Cycles/Year)") title("Periodogram of Relative Sunspot Number Data")
가산성 백색 잡음과 함께 각주파수 rad/sample과 rad/sample을 갖는 2개의 이산시간 정현파로 구성된 입력 신호의 주기도를 구합니다. rad/sample과 rad/sample에서 양측 주기도 추정값을 구합니다.
n = 0:319; x = cos(pi/4*n) + 0.5*sin(pi/2*n) + randn(size(n)); [pxx,w] = periodogram(x,[],[pi/4 pi/2]); pxx
pxx = 1×2
14.0589 2.8872
[pxx1,w1] = periodogram(x);
그 결과를 단측 주기도와 비교합니다. 얻어진 주기도 값은 단측 주기도 값의 1/2에 해당합니다. 특정 주파수 집합에서 주기도를 구하는 경우 출력값은 양측 추정값입니다.
plot(w1/pi,pxx1,w/pi,2*pxx,"o") legend("pxx1","2*pxx") xlabel("Normalized Frequency (\times \pi rad/sample)") title("Periodogram Power Spectral Density Estimate")
N(0,1) 가산성 백색 잡음에서 주파수 100Hz와 200Hz를 갖는 두 개의 사인파로 구성된 신호를 생성합니다. 샘플 레이트는 1kHz입니다. 100Hz와 200Hz에서 양측 주기도를 구합니다.
rng("default")
Fs = 1000;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t) + sin(2*pi*200*t) + randn(size(t));
freq = [100 200];
pxx = periodogram(x,[],freq,Fs)pxx = 1×2
0.2647 0.2313
이 예제에서는 주기도에서 신뢰한계를 사용하는 것을 다룹니다. 통계적 유의성에 있어 필요 조건은 아니지만, 주변 PSD 추정값의 신뢰 하한이 신뢰 상한을 초과하는 주기도의 주파수는 시계열에서 유의미한 진동을 명확하게 나타냅니다.
가산성 백색 N(0,1) 잡음에서 중첩된 형태의 100Hz 사인파와 150Hz 사인파로 구성된 신호를 생성합니다. 두 사인파의 진폭은 1입니다. 샘플 레이트는 1kHz입니다.
rng("default")
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + sin(2*pi*150*t) + randn(size(t));95% 신뢰한계를 갖는 주기도 PSD 추정값을 구합니다.
L = length(x); [pxx,f,pxxc] = periodogram(x,rectwin(L),L,Fs,ConfidenceLevel=0.95);
신뢰구간에 따라 주기도를 플로팅하고 100Hz와 150Hz 근처의 주파수 관심 영역을 확대합니다. 100Hz와 150Hz의 바로 인근에 있는 신뢰 하한은 100Hz와 150Hz 주변을 벗어난 신뢰 상한보다 상당히 높습니다.
plot(f,pow2db(pxx)) hold on plot(f,pow2db(pxxc),"-.",Color=[0.866 0.329 0]) xlim([85 175]) xlabel("Hz") ylabel("dB/Hz") title("Periodogram with 95%-Confidence Bounds")
가산성 잡음이 있는 100Hz 사인파의 주기도를 구합니다. 샘플 레이트는 1kHz입니다. "centered" 옵션을 사용하여 DC 기준 주기도를 구하고 결과를 플로팅합니다.
fs = 1000;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t)+randn(size(t));
periodogram(x,[],length(x),fs,"centered")
백색 가우스 잡음에 묻힌 200Hz 정현파로 구성된 신호를 생성합니다. 신호는 1초 동안 1kHz로 샘플링됩니다. 이 잡음의 분산은 0.01²입니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다.
rng("default")
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
N = length(t);
x = sin(2*pi*t*200) + 0.01*randn(size(t));FFT를 사용하여 신호 길이로 정규화된 신호의 파워 스펙트럼을 계산합니다. 정현파가 Bin 내부에 있으므로 모든 전력이 단일 주파수 샘플에 집중되어 있습니다. 단측 스펙트럼을 플로팅합니다. 피크의 인근을 확대합니다.
q = fft(x,N); ff = 0:Fs/N:Fs-Fs/N; ffts = q*q'/N^2
ffts = 0.4997
ff = ff(1:floor(N/2)+1);
q = q(1:floor(N/2)+1);
stem(ff,abs(q)/N,"*")
axis([190 210 0 0.55])
periodogram을 사용하여 신호의 파워 스펙트럼을 계산합니다. 핸(Hann) 윈도우를 지정하고 FFT 길이를 1024로 지정합니다. 200Hz에서 추정된 전력과 실제 값 사이의 백분율 차이를 구합니다.
wind = hann(N); [pun,fr] = periodogram(x,wind,1024,Fs,"power"); hold on stem(fr,pun)
periodogErr = abs(max(pun)-ffts)/ffts*100
periodogErr = 4.7349
파워 스펙트럼을 다시 계산하되, 이번에는 재할당을 사용합니다. 새 추정값을 플로팅하고 최댓값을 FFT 값과 비교합니다.
[pre,ft,pxx,fx] = periodogram(x,wind,1024,Fs,"power","reassigned"); stem(fx,pre) hold off legend(["Original" "Periodogram" "Reassigned"])
reassignErr = abs(max(pre)-ffts)/ffts*100
reassignErr = 0.0779
"power" 옵션을 사용하여 특정 주파수에서 정현파의 전력을 추정합니다.
1초 동안 1kHz로 샘플링된 100Hz 정현파를 생성합니다. 사인파의 진폭은 1.8이고, 이는 1.8²/2 = 1.62의 전력과 같습니다. "power" 옵션을 사용하여 전력을 추정합니다.
fs = 1000; t = 0:1/fs:1-1/fs; x = 1.8*cos(2*pi*100*t); [pxx,f] = periodogram(x,hamming(length(x)),length(x),fs,"power"); [pwrest,idx] = max(pxx); fprintf("The maximum power occurs at %3.1f Hz\n",f(idx))
The maximum power occurs at 100.0 Hz
fprintf("The power estimate is %2.2f\n",pwrest)The power estimate is 1.62
가산성 백색 가우스 잡음(AWGN)에서 3개의 정현파로 구성된 다중채널 신호의 샘플 1024개를 생성합니다. 정현파의 주파수는 rad/sample, rad/sample, rad/sample입니다. 주기도를 사용하여 신호의 PSD를 추정하고 플로팅합니다.
rng("default")
N = 1024;
n = 0:N-1;
w = pi./[2;3;4];
x = cos(w*n)' + randn(length(n),3);
periodogram(x)
윈도우를 사용하여 입력 신호의 수정된 주기도 PSD(파워 스펙트럼 밀도) 추정값을 반환하는 함수 myPeriodogram.m을 검토합니다. 함수는 이산 푸리에 변환 점의 개수를 입력 신호의 길이와 동일하게 지정합니다.
type myPeriodogramfunction [pxx,f] = myPeriodogram(inputData,window) %#codegen
nfft = length(inputData);
[pxx,f] = periodogram(inputData,window,nfft);
end
codegen (MATLAB Coder)을 사용하여 MEX 파일 을 생성합니다.
함수의
%#codegen지시문은 이 MATLAB® 코드가 코드 생성을 위한 코드임을 나타냅니다.-args옵션은 MEX 파일에 대한 입력값의 크기, 클래스, 실수/복소수 여부를 정의하는 샘플 인수를 지정합니다. 예를 들어,inputData를 1024×1 배정밀도 확률 벡터로 지정하고window를 길이가 1024인 해밍 윈도우로 지정합니다. 이후 MEX 함수를 호출할 때 1024개 샘플을 갖는 입력 신호와 윈도우를 사용하게 됩니다.MEX 함수에 다른 이름을 설정하려면
-o옵션을 사용하십시오.코드 생성 리포트를 생성하려면
codegen명령 끝에-report옵션을 추가하십시오.
codegen myPeriodogram -args {randn(1024,1),hamming(1024)}
Code generation successful.
periodogram 함수와 생성된 MEX 함수를 사용하여 1024개 샘플을 갖는 잡음이 있는 정현파의 PSD 추정값을 계산합니다. rad/sample의 정규화 정현파 주파수와 핸 윈도우를 지정합니다. 두 추정값을 플로팅하여 두 값이 일치하는 것을 확인합니다.
N = 1024; x = 2*cos(2*pi/5*(0:N-1)') + randn(N,1); periodogram(x,hann(N)) [pxMex,fMex] = myPeriodogram_mex(x,hann(N)); hold on plot(fMex/pi,pow2db(pxMex),"--") hold off grid on legend(["periodogram" "MEX function"])
R2026a 이후
4개의 신호에 대한 주기도 파워 스펙트럼 밀도(PSD) 추정값과 주기도 파워 스펙트럼을 지정된 대상 좌표축과 패널 컨테이너에 플로팅합니다.
3초 동안 10kHz의 샘플 레이트를 사용하여 4개의 진동 신호를 만듭니다.
Fs = 10e3;
t = 0:1/Fs:3;
x1 = sinc(Fs/2.5*(t-mean(t)));
x2 = sum(cos(2*pi*600*[1 3 5 7]'.*t),1) + randn(size(t))/1e4;
x3 = exp(1j*pi*sin(4*t)*Fs/10);
x4 = chirp(t,Fs/10,t(end),Fs/2.5,"quadratic");주기도 PSD 추정값과 파워 스펙트럼을 대상 좌표축에 플로팅하기
새 Figure 창의 남서쪽 코너와 북동쪽 코너에 두 개의 좌표축을 만듭니다.
fig = figure; ax1 = axes(fig,Position=[0.09 0.1 0.60 0.45]); ax2 = axes(fig,Position=[0.50 0.7 0.47 0.25]);
신호 x1과 x2의 주기도 PSD 추정값과 파워 스펙트럼을 각각 Figure의 남서쪽 좌표축과 북동쪽 좌표축에 플로팅합니다. 카이저 윈도우와 512개 DFT 점을 사용합니다.
g = kaiser(length(t),5);
nfft = 512;
periodogram(x1,g,nfft,Fs,Parent=ax1)
periodogram(x2,g,nfft,Fs,"power",Parent=ax2)
주기도 PSD 추정값을 대상 UI 좌표축에 플로팅하기
새 UI Figure 창의 북서쪽 코너에 좌표축을 만듭니다.
uif = uifigure(Position=[100 100 720 540]); ax3 = uiaxes(uif,Position=[5 305 300 200]);
신호 x3에 대한 주기도 PSD 추정값을 Figure 좌표축에 플로팅합니다. 중심이 0kHz에 맞춰진 주파수를 표시합니다.
periodogram(x3,g,nfft,Fs,"centered",Parent=ax3) title(ax3,"Periodogram in UI Axes")
재할당된 주기도 PSD 추정값을 대상 패널 컨테이너에 플로팅하기
UI Figure 창의 남동쪽 코너에 패널 컨테이너를 추가합니다.
p = uipanel(uif,Position=[300 5 400 325], ... Title="Periodogram in Panel Container", ... BackgroundColor="white");
신호 x4에 대한 재할당된 주기도 PSD 추정값을 패널 컨테이너에 플로팅합니다.
periodogram(x4,g,nfft,Fs,"reassigned",Parent=p)
입력 인수
출력 인수
세부 정보
참고 문헌
[1] Auger, François, and Patrick Flandrin. "Improving the Readability of Time-Frequency and Time-Scale Representations by the Reassignment Method." IEEE® Transactions on Signal Processing. Vol. 43, May 1995, pp. 1068–1089.
[2] Fulop, Sean A., and Kelly Fitz. "Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications." Journal of the Acoustical Society of America. Vol. 119, January 2006, pp. 360–371.
