pwelch
Welch의 파워 스펙트럼 밀도 추정값
구문
설명
[는 신호 pxx,f] = pwelch(x,win,nOverlap,freqSpec)x에 대한 Welch PSD 추정값과 주파수 f(rad/sample)를 반환합니다. 여기서 pwelch 함수는 다음을 수행합니다.
win을 사용하여 신호를 여러 세그먼트로 나누고 윈도우를 적용합니다.nOverlap에 지정된 중첩 길이를 사용하여 인접 세그먼트 간에 샘플을 중첩합니다.freqSpec에 지정된 DFT 지점 개수 또는 주파수에서, 윈도우가 적용된 각 세그먼트에 대해 이산 푸리에 변환(DFT)을 계산합니다.
이러한 입력 인수에 디폴트 값을 사용하려면 빈 값 []을 지정하십시오.
[___] = pwelch(___,는 위에 열거된 구문에 주파수 범위, 스펙트럼 유형, 추적 모드도 지정합니다. 이러한 입력 인수의 조합을 지정할 수 있습니다.freqRange,spectrumType,trace)
pwelch(___)를 출력 인수를 지정하지 않고 사용하면 Welch PSD 추정값 또는 파워 스펙트럼을 현재 Figure 창에 플로팅합니다.
예제
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 Welch PSD 추정값을 구합니다.
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 사인파를 생성합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 신호의 길이는 개 샘플입니다.
rng("default")
n = 0:319;
x = cos(pi/4*n) + randn(size(n));디폴트 해밍 윈도우와 DFT 길이를 사용하여 Welch PSD 추정값을 구합니다. 디폴트 세그먼트 길이는 71개 샘플이고 DFT 길이는 256개 점이며, 이 경우 rad/sample의 주파수 분해능이 생성됩니다. 신호가 실수 값이므로 주기도는 단측이고 256/2+1개 점이 있습니다. Welch PSD 추정값을 플로팅합니다.
pxx = pwelch(x); pwelch(x)
계산을 반복합니다.
신호를 길이가 인 섹션으로 나눕니다. 이 작업은 50% 중첩을 갖는 세그먼트의 개수가 8에 최대한 근접하되 초과하지는 않도록 하는 가능한 한 가장 긴 섹션으로 나누는 것과 동일합니다.
해밍 윈도우를 사용하여 섹션에 윈도우를 적용합니다.
연속 섹션 간에 50% 중첩을 지정합니다.
개의 지점을 사용하여 FFT를 계산합니다. 여기서 입니다.
두 접근 방식이 동일한 결과를 제공하는지 확인합니다.
Nx = length(x); nsc = floor(Nx/4.5); nov = floor(nsc/2); nff = max(256,2^nextpow2(nsc)); pxxt = pwelch(x,hamming(nsc),nov,nff); maxerr = max(abs(abs(pxxt(:)) - abs(pxx(:))))
maxerr = 0
섹션 간에 50% 중첩을 가지면서 길이가 동일한 8개 섹션으로 신호를 나눕니다. 이전 단계에서와 동일한 FFT 길이를 지정합니다. 단시간 Welch PSD 추정값을 계산하고 이전 두 절차와 같은 결과를 제공하는지 확인합니다.
ns = 8; ov = 0.5; lsc = floor(Nx/(ns-(ns-1)*ov)); pxxt8 = pwelch(x,lsc,floor(ov*lsc),nff); maxerr8 = max(abs(abs(pxxt8(:)) - abs(pxx(:))))
maxerr8 = 0
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 Welch PSD 추정값을 구합니다.
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 사인파를 생성합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 신호 샘플은 512개입니다.
rng("default")
n = 0:511;
x = cos(pi/3*n) + randn(size(n));신호를 132개 샘플 길이의 세그먼트로 나누어 Welch PSD 추정값을 구합니다. 신호 세그먼트에 132개 샘플 길이의 해밍 윈도우가 곱해집니다. 중첩된 샘플의 개수가 지정되지 않아 132/2 = 66으로 설정됩니다. DFT 길이는 256개 점입니다. 이 경우 rad/sample의 주파수 분해능이 생성됩니다. 신호가 실수 값이므로 PSD 추정값은 단측이고 256/2+1 = 129개 점이 있습니다. PSD를 정규화 주파수의 함수로 플로팅합니다.
segmentLength = 132; [pxx,w] = pwelch(x,segmentLength); plot(w/pi,pow2db(pxx)) xlabel("Normalized Frequency (\times \pi rad/sample)") title("Welch Power Spectral Density Estimate")
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 Welch PSD 추정값을 구합니다.
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 사인파를 생성합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 신호의 길이는 320개 샘플입니다.
rng("default")
n = 0:319;
x = cos(pi/4*n) + randn(size(n));신호를 100개 샘플 길이의 세그먼트로 나누어 Welch PSD 추정값을 구합니다. 신호 세그먼트에 100개 샘플 길이의 해밍 윈도우가 곱해집니다. 중첩된 샘플의 개수는 25입니다. DFT 길이는 256개 점입니다. 이 경우 rad/sample의 주파수 분해능이 생성됩니다. 신호가 실수 값이므로 PSD 추정값은 단측이고 256/2+1개 점이 있습니다.
segmentLength = 100; noverlap = 25; pxx = pwelch(x,segmentLength,noverlap); plot(pow2db(pxx)) xlabel("Frequency Samples") title("Welch Power Spectral Density Estimate")
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 Welch PSD 추정값을 구합니다.
가산성 백색 잡음과 함께 각주파수 rad/sample을 갖는 사인파를 생성합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 신호의 길이는 320개 샘플입니다.
rng("default")
n = 0:319;
x = cos(pi/4*n) + randn(size(n));신호를 100개 샘플 길이의 세그먼트로 나누어 Welch PSD 추정값을 구합니다. 디폴트 중첩 50%를 사용합니다. 주파수 rad/sample이 DFT Bin(Bin 81)에 대응되도록 DFT 길이를 640개 점으로 지정합니다. 신호가 실수 값이므로 PSD 추정값은 단측이고 640/2+1개 점이 있습니다.
segmentLength = 100; nfft = 640; pxx = pwelch(x,segmentLength,[],nfft); plot(pow2db(pxx)) xlabel("Frequency Points") title("Welch Power Spectral Density Estimate")
가산성 N(0,1) 백색 잡음에서 100Hz 정현파로 구성된 신호를 생성합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 샘플 레이트는 1kHz이고 신호의 지속 시간은 5초입니다.
rng("default")
Fs = 1000;
t = 0:1/Fs:5-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));위에 나온 신호의 Welch의 중첩 세그먼트 평균화 PSD 추정값을 구합니다. 세그먼트 길이는 300개의 중첩된 샘플을 가지는 500개 샘플을 사용합니다. 100Hz가 정확히 DFT Bin에 떨어지도록 500개 DFT 점을 사용합니다. Hz 단위의 주파수 벡터를 출력하도록 샘플 레이트를 입력합니다. 결과를 플로팅합니다.
[pxx,f] = pwelch(x,500,300,500,Fs); plot(f,pow2db(pxx)) xlabel("Frequency (Hz)") ylabel("PSD (dB/Hz)")
0.1초 동안 200kHz로 샘플링된, 3개의 잡음이 있는 정현파와 하나의 처프(Chirp)로 구성된 신호를 생성합니다. 정현파의 주파수는 1kHz, 10kHz, 20kHz입니다. 정현파는 각기 다른 진폭과 잡음 수준을 가집니다. 잡음이 없는 처프는 샘플링 동안 20kHz에서 시작하여 30kHz까지 선형적으로 증가하는 주파수를 가집니다.
rng("default")
Fs = 200e3;
Fc = [1 10 20]'*1e3;
Ns = 0.1*Fs;
t = (0:Ns-1)/Fs;
x = [1 1/10 10]*sin(2*pi*Fc*t) + [1/200 1/2000 1/20]*randn(3,Ns);
x = x + chirp(t,20e3,t(end),30e3);Welch PSD 추정값과 신호에 대한 최댓값 유지 스펙트럼 및 최솟값 유지 스펙트럼을 계산합니다. 결과를 플로팅합니다.
[pxx,f] = pwelch(x,[],[],[],Fs); pmax = pwelch(x,[],[],[],Fs,"maxhold"); pmin = pwelch(x,[],[],[],Fs,"minhold"); plot(f,pow2db(pxx)) hold on plot(f,pow2db([pmax pmin]),":") hold off xlabel("Frequency (Hz)") ylabel("PSD (dB/Hz)") legend(["pwelch" "maxhold" "minhold"])
절차를 반복하되, 이번에는 중심에 있는 파워 스펙트럼 추정값을 계산합니다.
[pxx,f] = pwelch(x,[],[],[],Fs,"centered","power"); pmax = pwelch(x,[],[],[],Fs,"maxhold","centered","power"); pmin = pwelch(x,[],[],[],Fs,"minhold","centered","power"); plot(f,pow2db(pxx)) hold on plot(f,pow2db([pmax pmin]),":") hold off xlabel("Frequency (Hz)") ylabel("Power (dB)") legend(["pwelch" "maxhold" "minhold"])
이 예제에서는 Welch의 중첩 세그먼트 평균화(WOSA) PSD 추정값과 함께 신뢰한계를 사용하는 경우를 다룹니다. 통계적 유의성에 있어 필요 조건은 아니지만, 주변 PSD 추정값의 신뢰 하한이 신뢰 상한을 초과하는 Welch 추정값의 주파수는 시계열에서 유의미한 진동을 명확하게 나타냅니다.
가산성 백색 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% 신뢰한계를 갖는 WOSA 추정값을 구합니다. 세그먼트 길이를 200으로 설정하고 세그먼트를 50%(100개 샘플) 중첩합니다.
L = 200; win = hamming(L); nOverlap = 100; [pxx,f,pxxc] = pwelch(x,win,nOverlap,200,Fs,ConfidenceLevel=0.95);
신뢰구간에 따라 WOSA PSD 추정값을 플로팅하고 100Hz와 150Hz 근처의 주파수 관심 영역을 확대합니다. 100Hz와 150Hz의 바로 인근에 있는 신뢰 하한은 100Hz와 150Hz 주변을 벗어난 신뢰 상한보다 상당히 높습니다.
plot(f,pow2db(pxx)) hold on plot(f,pow2db(pxxc),"-.",Color=[0.866 0.329 0]) hold off xlim([25 250]) xlabel("Frequency (Hz)") ylabel("PSD (dB/Hz)") title("Welch Estimate with 95%-Confidence Bounds")
가산성 백색 잡음에서 100Hz 정현파로 구성된 신호를 생성합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 샘플 레이트는 1kHz이고 신호의 지속 시간은 5초입니다.
rng("default")
Fs = 1000;
t = 0:1/Fs:5-1/Fs;
noisevar = 1/4;
x = cos(2*pi*100*t) + sqrt(noisevar)*randn(size(t));Welch 방법을 사용하여 DC 기준 파워 스펙트럼을 구합니다. 세그먼트 길이는 300개의 중첩된 샘플을 가지는 500개 샘플을, DFT 길이는 500개 점을 사용합니다.
[pxx,f] = pwelch(x,500,300,500,Fs,"centered","power");
결과를 플로팅합니다. -100Hz와 100Hz에서의 전력은 진폭이 1인 실수 값 사인파에 대해 예상되는 전력 1/4에 가깝습니다. 1/4에서의 편차는 가산성 잡음의 효과로 인한 것입니다.
plot(f,pow2db(pxx)) xlabel("Frequency (Hz)") ylabel("Magnitude (dB)") grid
가산성 백색 가우스 잡음(AWGN)에서 3개의 정현파로 구성된 다중채널 신호의 샘플 1024개를 생성합니다. 정현파의 주파수는 rad/sample, rad/sample, rad/sample입니다. Welch 방법을 사용하여 신호의 PSD를 추정하고 플로팅합니다.
N = 1024;
n = 0:N-1;
w = pi./[2;3;4];
rng("default")
x = cos(w*n)' + randn(length(n),3);
pwelch(x)
R2026a 이후
4개의 신호에 대한 Welch 파워 스펙트럼 밀도(PSD) 추정값과 Welch 파워 스펙트럼을 지정된 대상 좌표축과 패널 컨테이너에 플로팅합니다.
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");Welch PSD 추정값과 파워 스펙트럼을 대상 좌표축에 플로팅하기
새 Figure 창의 남서쪽 코너와 북동쪽 코너에 두 개의 좌표축을 만듭니다.
fig = figure; ax1 = axes(fig,Position=[0.09 0.1 0.52 0.45]); ax2 = axes(fig,Position=[0.55 0.7 0.42 0.25]);
신호 x1과 x2의 Welch PSD 추정값과 Welch 파워 스펙트럼을 각각 Figure의 남서쪽 좌표축과 북동쪽 좌표축에 플로팅합니다. 256개 샘플의 카이저 윈도우, 220개 샘플의 중첩 길이, 512개 DFT 점을 사용합니다.
g = kaiser(256,5);
ol = 220;
nfft = 512;
pwelch(x1,g,ol,nfft,Fs,Parent=ax1)
pwelch(x2,g,ol,nfft,Fs,"power",Parent=ax2)
Welch PSD 추정값을 대상 UI 좌표축에 플로팅하기
새 UI Figure 창의 북서쪽 코너에 좌표축을 만듭니다.
uif = uifigure(Position=[100 100 720 540]); ax3 = uiaxes(uif,Position=[5 305 300 200]);
신호 x3에 대한 Welch PSD 추정값을 Figure 좌표축에 플로팅합니다. 중심이 0kHz에 맞춰진 주파수를 표시합니다.
pwelch(x3,g,ol,nfft,Fs,"centered",Parent=ax3) title(ax3,"Welch PSD in UI Axes")
Welch PSD 추정값을 대상 패널 컨테이너에 플로팅하기
UI Figure 창의 남동쪽 코너에 패널 컨테이너를 추가합니다.
p = uipanel(uif,Position=[300 5 400 325], ... Title="Welch PSD in Panel Container", ... BackgroundColor="white");
신호 x4에 대한 Welch PSD 추정값을 패널 컨테이너에 플로팅합니다. 신뢰수준을 90%로 설정합니다.
pwelch(x4,g,ol,nfft,Fs,ConfidenceLevel=0.9,Parent=p)
입력 인수
출력 인수
세부 정보
참고 문헌
[1] Hayes, Monson H. Statistical Digital Signal Processing and Modeling. New York: John Wiley & Sons, 1996.
[2] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.
