pmtm

멀티테이퍼 파워 스펙트럼 밀도 추정값

구문

pxx = pmtm(x)

pxx = pmtm(x,'Tapers',tapertype)

pxx = pmtm(x,nw)

pxx = pmtm(x,m,'Tapers','sine')

pxx = pmtm(___,nfft)

[pxx,w] = pmtm(___)

[pxx,f] = pmtm(___,fs)

[pxx,w] = pmtm(x,nw,w)

[pxx,w] = pmtm(x,m,'Tapers','sine',w)

[pxx,f] = pmtm(___,f,fs)

[___] = pmtm(___,freqrange)

[___,pxxc] = pmtm(___,'ConfidenceLevel',probability)

[___] = pmtm(___,'DropLastTaper',dropflag)

[___] = pmtm(___,method)

[___] = pmtm(x,e,v,___)

[___] = pmtm(x,dpss_params,___)

pmtm(___)

설명

pxx = pmtm(x)는 이산 장형 타원체(슬레피안) 시퀀스를 테이퍼로 사용하여 입력 신호 x의 Thomson 멀티테이퍼 파워 스펙트럼 밀도(PSD) 추정값 pxx를 반환합니다.

pxx = pmtm(x,'Tapers',tapertype)은 멀티테이퍼 PSD 추정값을 계산할 때 사용할 테이퍼의 유형을 지정합니다. 'Tapers',tapertype 이름-값 쌍은 함수 호출에서 x 이후의 아무 위치에나 지정할 수 있습니다.

pxx = pmtm(x,nw)는 슬레피안 테이퍼를 사용하여 PSD 추정값을 계산할 때 시간과 반대역폭의 곱 nw를 사용하여 주파수 분해능을 제어합니다.

pxx = pmtm(x,m,'Tapers','sine')은 사인 테이퍼를 사용하여 PSD 추정값을 계산할 때 적용할 테이퍼의 개수 또는 평균 가중치를 지정합니다.

pxx = pmtm(___,nfft)는 위에 열거된 구문 중 하나와 함께 nfft개의 이산 푸리에 변환(DFT) 점을 사용합니다. nfft가 신호 길이보다 크면 x가 길이 nfft까지 0으로 채워집니다. nfft가 신호 길이보다 작으면 신호가 모듈로 nfft로 순환됩니다.

[pxx,w] = pmtm(___)은 pxx가 계산되는 지점의 정규화 주파수로 구성된 벡터를 반환합니다.

[pxx,f] = pmtm(___,fs)는 주파수 벡터 f를 반환합니다(단위: 단위 시간당 주기). fs는 함수 호출에서 x, nw(사인 테이퍼의 경우 m), nfft 뒤에 와야 합니다. 샘플 레이트를 입력하고 위에 열거된 인수의 디폴트 값을 그대로 사용하려면 이 인수를 빈 값 []로 지정하십시오.

[pxx,w] = pmtm(x,nw,w)는 w에 지정된 정규화 주파수에서 슬레피안 시퀀스를 사용하여 계산된 멀티테이퍼 PSD 추정값을 반환합니다. 벡터 w는 적어도 2개의 요소를 가져야 합니다.

[pxx,w] = pmtm(x,m,'Tapers','sine',w)는 w에 지정된 정규화 주파수에서 사인 테이퍼를 사용하여 계산된 멀티테이퍼 PSD 추정값을 반환합니다. 벡터 w는 적어도 2개의 요소를 가져야 합니다.

[pxx,f] = pmtm(___,f,fs)는 f에 지정된 주파수에서 멀티테이퍼 PSD 추정값을 계산합니다. 벡터 f는 샘플 레이트 fs와 동일한 단위의 요소를 적어도 2개 가져야 합니다.

[___] = pmtm(___,freqrange)는 freqrange로 지정된 주파수 범위에 대한 멀티테이퍼 PSD 추정값을 반환합니다.

[___,pxxc] = pmtm(___,'ConfidenceLevel',probability)는 PSD 추정값의 probability × 100% 신뢰구간을 pxxc로 반환합니다.

[___] = pmtm(___,'DropLastTaper',dropflag)는 멀티테이퍼 PSD 추정값을 계산할 때 pmtm이 마지막 슬레피안 테이퍼를 제외하는지 여부를 지정합니다.

[___] = pmtm(___,method)는 테이퍼가 적용된 개별 PSD 추정값을 method에 지정된 메서드를 사용하여 결합합니다. 이 구문은 슬레피안 테이퍼에만 적용됩니다.

[___] = pmtm(x,e,v,___)는 e의 슬레피안 테이퍼와 v의 고유값을 사용하여 PSD를 계산합니다. dpss를 사용하여 e와 v를 구합니다.

[___] = pmtm(x,dpss_params,___)는 셀형 배열 dpss_params를 사용하여 입력 인수를 dpss에 전달합니다. 이 구문은 슬레피안 테이퍼에만 적용됩니다.

pmtm(___)에 출력 인수를 지정하지 않으면 현재 Figure 창에 멀티테이퍼 PSD 추정값을 플로팅합니다.

예제

디폴트 입력값을 사용한 멀티테이퍼 추정값

가산성 N(0,1) 백색 잡음과 함께 각주파수 $$\pi /4$$ rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 멀티테이퍼 PSD 추정값을 구합니다.

가산성 N(0,1) 백색 잡음과 함께 각주파수 $$\pi /4$$ rad/sample을 갖는 사인파를 생성합니다. 신호의 길이는 320개 샘플입니다. 시간과 반대역폭의 곱의 디폴트 값 4와 DFT 길이를 사용하여 멀티테이퍼 PSD 추정값을 구합니다. DFT 점의 디폴트 개수는 512개입니다. 신호가 실수 값이므로 PSD 추정값은 단측이고 PSD 추정값에 512/2+1개 점이 있습니다.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
pxx = pmtm(x);

멀티테이퍼 PSD 추정값을 플로팅합니다.

pmtm(x)

사인 테이퍼와 슬레피안 테이퍼

가산성 N(0,1) 백색 가우스 잡음에 묻힌 2채널 신호의 샘플 2048개를 생성합니다.

• 첫 번째 채널은 정규화 주파수가 π/3 rad/sample 및 π/5 rad/sample인 두 개의 정현파로 구성되어 있습니다. 첫 번째 정현파의 진폭은 두 번째 정현파의 두 배입니다.

• 두 번째 채널의 정규화 주파수는 π/4 rad/sample입니다.

멀티테이퍼 방법을 사용하여 0.1π rad/sample에서 0.4π rad/sample까지 1024개 샘플 구간에 대해 신호의 PSD를 추정합니다. 동일한 가중치가 적용된 13개의 사인 테이퍼를 사용합니다.

n = (0:2047)';

x = [sin(pi./[3 5].*n)*[2 1]' sin(pi/4*n)] + randn(length(n),2);

w = linspace(0.1,0.4,1024);

ntp = 13;
pmtm(x,ntp,'Tapers','sine',w*pi)

이번에는 13개의 테이퍼에 선형 내림차순으로 가중치를 적용하여 계산을 반복합니다. 'Tapers','sine' 이름-값 쌍은 함수 호출에서 x 이후의 아무 위치에나 지정할 수 있습니다.

pmtm(x,(ntp:-1:1)/sum(1:ntp),w*pi,'Tapers','sine')

이번에는 13개의 슬레피안 테이퍼를 사용하고 시간과 반대역폭의 곱을 7.5로 지정하여 계산을 반복합니다.

nw = 7.5;

pmtm(x,{nw,ntp},w*pi)

이번에는 샘플 레이트를 2kHz로 지정하여 계산을 반복합니다.

fs = 2e3;

pmtm(x,{nw,ntp},w*(fs/2),fs)

시간과 반대역폭의 곱 지정하기

지정된 시간과 반대역폭의 곱으로 멀티테이퍼 PSD 추정값을 구합니다.

가산성 N(0,1) 백색 잡음과 함께 각주파수 $$\pi /4$$ rad/sample을 갖는 사인파를 생성합니다. 신호의 길이는 320개 샘플입니다. 시간과 반대역폭의 곱을 2.5로 지정하여 멀티테이퍼 PSD 추정값을 구합니다. 분해능 대역폭은 $$[-2.5\pi /320,2.5\pi /320]$$ rad/sample입니다. DFT 점의 디폴트 개수는 512개입니다. 신호가 실수 값이므로 PSD 추정값은 단측이고 PSD 추정값에 512/2+1개 점이 있습니다.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
pmtm(x,2.5)

신호 길이와 동일한 DFT 길이

가산성 N(0,1) 백색 잡음과 함께 각주파수 $$\pi /4$$ rad/sample을 갖는 이산시간 정현파로 구성된 입력 신호의 멀티테이퍼 PSD 추정값을 구합니다. 신호 길이와 동일한 DFT 길이를 사용합니다.

가산성 N(0,1) 백색 잡음과 함께 각주파수 $$\pi /4$$ rad/sample을 갖는 사인파를 생성합니다. 신호의 길이는 320개 샘플입니다. 시간과 반대역폭의 곱을 3으로 지정하고 DFT 길이를 신호 길이와 동일하게 하여 멀티테이퍼 PSD 추정값을 구합니다. 신호가 실수 값을 가지기 때문에 기본적으로 길이가 320/2+1인 단측 PSD 추정값이 반환됩니다.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
pmtm(x,3,length(x))

샘플 레이트를 사용한 멀티테이퍼 추정값

1kHz로 샘플링된 신호의 멀티테이퍼 PSD 추정값을 구합니다. 이 신호는 가산성 N(0,1) 백색 잡음이 있는 100Hz의 사인파입니다. 신호의 지속 시간은 2s입니다. 시간과 반대역폭의 곱을 3으로 지정하고 DFT 길이를 신호 길이와 동일하게 합니다.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));
[pxx,f] = pmtm(x,3,length(x),fs);

멀티테이퍼 PSD 추정값을 플로팅합니다.

pmtm(x,3,length(x),fs)

단위 가중치를 사용한 평균 싱글테이퍼 추정값

테이퍼가 적용된 직접 스펙트럼 추정값 각각에 평균적으로 동일한 가중치를 적용하여, 멀티테이퍼 PSD 추정값을 구합니다.

1kHz로 샘플링된 신호의 멀티테이퍼 PSD 추정값을 구합니다. 이 신호는 가산성 N(0,1) 백색 잡음이 있는 100Hz의 사인파입니다. 신호의 지속 시간은 2s입니다. 시간과 반대역폭의 곱을 3으로 지정하고 DFT 길이를 신호 길이와 동일하게 합니다. 'unity' 옵션을 사용하여, 테이퍼가 적용된 직접 스펙트럼 추정값 각각에 평균적으로 동일한 가중치를 적용합니다.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));
[pxx,f] = pmtm(x,3,length(x),fs,'unity');

멀티테이퍼 PSD 추정값을 플로팅합니다.

pmtm(x,3,length(x),fs,'unity')

DPSS 시퀀스와 이에 대응하는 주파수 영역 집중도

이 예제에서는 DPSS 시퀀스의 주파수 영역 집중도를 살펴봅니다. 또한 슬레피안 시퀀스를 미리 계산하고 그중에서 분해능 대역폭에 집중된 에너지가 99%를 넘는 시퀀스만 선택하여, 입력 신호의 멀티테이퍼 PSD 추정값을 구합니다.

이 신호는 가산성 N(0,1) 백색 잡음이 있는 100Hz의 사인파입니다. 신호의 지속 시간은 2초입니다.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));

시간과 반대역폭의 곱을 3.5로 설정합니다. 신호 길이가 2000개 샘플이고 샘플링 간격이 0.001초인 경우 분해능 대역폭은 [-1.75,1.75]Hz가 됩니다. 처음 10개의 슬레피안 시퀀스를 계산하고, 지정된 분해능 대역폭에서 해당 주파수 집중도를 살펴봅니다.

[e,v] = dpss(length(x),3.5,10);
lv = length(v);

stem(1:lv,v,"filled")
ylim([0 1.2])
title("Proportion of Energy in [-w,w] of k-th Slepian Sequence")

에너지 집중도가 99%가 넘는 슬레피안 시퀀스 개수를 확인합니다. 선택한 DPSS 시퀀스를 사용하여 멀티테이퍼 PSD 추정값을 구합니다. 선택한 테이퍼를 모두 사용하려면 DropLastTaper를 false로 설정하십시오.

hold on
plot([1 lv],0.99*[1 1])
hold off

idx = find(v>0.99,1,'last')

idx = 
5

[pxx,f] = pmtm(x,e(:,1:idx),v(1:idx),length(x),fs,DropLastTaper=false);

멀티테이퍼 PSD 추정값을 플로팅합니다.

pmtm(x,e(:,1:idx),v(1:idx),length(x),fs,DropLastTaper=false)

DC 기준 멀티테이퍼 PSD 추정값

가산성 N(0,1) 잡음이 있는 100Hz 사인파의 멀티테이퍼 PSD 추정값을 구합니다. 데이터는 1kHz로 샘플링됩니다. 'centered' 옵션을 사용하여 DC 기준 PSD를 구합니다.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));
[pxx,f] = pmtm(x,3.5,length(x),fs,'centered');

DC 기준 PSD 추정값을 플로팅합니다.

pmtm(x,3.5,length(x),fs,'centered')

95% 신뢰한계의 상한값과 하한값

다음 예제에서는 멀티테이퍼 PSD 추정값에서 신뢰한계를 사용하는 것을 다룹니다. 통계적 유의성에 있어 필요 조건은 아니지만, 주변 PSD 추정값의 신뢰 하한이 신뢰 상한을 초과하는 멀티테이퍼 PSD 추정값의 주파수는 시계열에서 유의미한 진동을 명확하게 나타냅니다.

가산성 백색 N(0,1) 잡음에서 중첩된 형태의 100Hz 사인파와 150Hz 사인파로 구성된 신호를 생성합니다. 두 사인파의 진폭은 1입니다. 샘플링 주파수는 1kHz입니다. 신호의 지속 시간은 2s입니다.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+cos(2*pi*150*t)+randn(size(t));

95% 신뢰한계를 갖는 멀티테이퍼 PSD 추정값을 구합니다. 신뢰구간에 따라 PSD 추정값을 플로팅하고 100Hz와 150Hz 근처의 주파수 관심 영역을 확대합니다.

[pxx,f,pxxc] = pmtm(x,3.5,length(x),fs,'ConfidenceLevel',0.95);

plot(f,10*log10(pxx))
hold on
plot(f,10*log10(pxxc),'r-.')
xlim([85 175])
xlabel('Hz')
ylabel('dB')
title('Multitaper PSD Estimate with 95%-Confidence Bounds')

100Hz와 150Hz의 바로 인근에 있는 신뢰 하한은 100Hz와 150Hz 주변을 벗어난 신뢰 상한보다 상당히 높습니다.

다중채널 신호의 멀티테이퍼 PSD 추정값

가산성 $$N(0,1)$$ 백색 가우스 잡음(AWGN)에서 3개의 정현파로 구성된 다중채널 신호의 샘플 1024개를 생성합니다. 정현파의 주파수는 $$\pi /2$$ rad/sample, $$\pi /3$$ rad/sample, $$\pi /4$$ rad/sample입니다. Thomson의 멀티테이퍼 방법을 사용하여 신호의 PSD를 추정하고 플로팅합니다.

N = 1024;
n = 0:N-1;

w = pi./[2;3;4];
x = cos(w*n)' + randn(length(n),3);

pmtm(x)

입력 인수

x — 입력 신호

벡터 | 행렬

입력 신호로, 행 벡터나 열 벡터 또는 행렬로 지정됩니다. x가 행렬이면, 이 행렬의 각 열은 독립적인 채널로 처리됩니다.

예: cos(pi/4*(0:159))+randn(1,160)은 단일채널 행 벡터 신호입니다.

예: cos(pi./[4;2]*(0:159))'+randn(160,2)는 2채널 신호입니다.

데이터형: single | double
복소수 지원 여부: 예

tapertype — 테이퍼 유형

'slepian' (디폴트 값) | 'sine'

테이퍼 유형으로, 'slepian' 또는 'sine'으로 지정됩니다.

• 'slepian' — 테이퍼로 이산 장형 타원체(슬레피안) 시퀀스를 사용합니다.

• 'sine' — 사인 테이퍼를 사용합니다.

'Tapers',tapertype 이름-값 쌍은 함수 호출에서 x 이후의 아무 위치에나 지정할 수 있습니다.

데이터형: char | string

nw — 시간과 반대역폭의 곱

4 (디폴트 값) | 양의 스칼라

시간과 반대역폭의 곱으로, 양의 스칼라로 지정됩니다. pmtm 함수는 PSD 추정값에 2 × nw – 1개 슬레피안 테이퍼를 사용합니다. nw로는 일반적으로 2, 5/2, 3, 7/2 또는 4를 사용합니다.

멀티테이퍼 스펙트럼 추정 시 사용자는 멀티테이퍼 추정값의 분해능 대역폭 [–W,W]를 지정합니다. 여기서 어느 정도 작은 k > 1에 대해 W = k/NΔt입니다. 이와 동등한 의미에서, W는 DFT 주파수 분해능의 어느 정도 작은 배수입니다. 시간과 반대역폭의 곱은 분해능 반대역폭과 입력 신호에 포함된 샘플 개수 N을 곱한 것입니다. 푸리에 변환이 [–W,W]에 잘 집중되어 있는(고유값이 단위 값에 가까움) 슬레피안 테이퍼 개수는 2NW – 1입니다.

m — 사인 테이퍼 개수 또는 평균 가중치

7 (디폴트 값) | 정수 스칼라 | 벡터

사인 테이퍼 개수 또는 평균 가중치로, 정수 스칼라 또는 벡터로 지정됩니다.

• m이 스칼라인 경우 PSD 추정값을 계산할 때 데이터 윈도우로 사용되는 사인 테이퍼의 개수를 나타냅니다. 사인 테이퍼에는 가중치가 균일하게 적용됩니다.

• m이 벡터인 경우 PSD 추정값을 계산할 때 사인 테이퍼의 평균을 계산하는 데 사용되는 가중치를 나타냅니다. m의 길이는 사용할 테이퍼의 개수를 나타냅니다. m의 요소는 합이 1이어야 합니다.

데이터형: single | double

nfft — DFT를 적용할 지점의 개수

max(256,2^nextpow2(length(x))) (디폴트 값) | 정수 | []

DFT를 적용할 지점의 개수로, 양의 정수로 지정됩니다. 실수 값을 갖는 입력 신호 x의 경우, nfft가 짝수이면 PSD 추정값 pxx의 길이가 (nfft/2 + 1)이고, nfft가 홀수이면 길이가 (nfft + 1)/2입니다. 복소수 값을 갖는 입력 신호 x의 경우, PSD 추정값은 항상 길이가 nfft입니다. nfft가 빈 값으로 지정되면 디폴트 nfft가 사용됩니다.

데이터형: single | double

fs — 샘플 레이트

양의 스칼라

샘플 레이트로, 양의 스칼라로 지정됩니다. 샘플 레이트는 단위 시간당 샘플 개수입니다. 시간 단위가 초이면 샘플 레이트의 단위는 Hz입니다.

w — 정규화 주파수

벡터

정규화 주파수로, 적어도 2개의 요소를 가진 행 벡터나 열 벡터로 지정됩니다. 정규화 주파수는 rad/sample을 단위로 사용합니다.

예: w = [pi/4 pi/2]

데이터형: double | single

f — 주파수

벡터

주파수로, 적어도 2개의 요소를 가진 행 벡터나 열 벡터로 지정됩니다. 주파수는 단위 시간당 주기를 단위로 사용합니다. 단위 시간은 샘플 레이트 fs로 지정됩니다. fs의 단위가 샘플/초이면 f는 Hz 단위입니다.

예: fs = 1000; f = [100 200]

데이터형: double | single

dropflag — 마지막 DPSS 시퀀스를 제외할지 또는 유지할지 여부를 나타내는 플래그

true (디폴트 값) | false

마지막 DPSS 시퀀스를 제외할지 또는 유지할지 여부를 나타내는 플래그로, 논리형으로 지정됩니다. 디폴트 값은 true이며 pmtm은 마지막 테이퍼를 제외시킵니다. 멀티테이퍼 추정값에서 처음 2NW – 1개 DPSS 시퀀스는 단위 값에 가까운 고유값을 갖습니다. 2NW – 1개 보다 작은 시퀀스를 사용하는 경우, 모든 테이퍼의 고유값이 1에 가까울 수 있으며 마지막 테이퍼를 유지하려면 dropflag를 false로 지정하면 됩니다.

method — 테이퍼가 적용된 개별 PSD 추정값에 대한 가중치

'adapt' (디폴트 값) | 'eigen' | 'unity'

테이퍼가 적용된 개별 PSD 추정값에 대한 가중치로, 'adapt', 'eigen' 또는 'unity' 중 하나로 지정됩니다. 디폴트 값은 Thomson의 주파수 종속적 가변 가중치인 'adapt'입니다. 이러한 가중치 계산 방법은 [2]의 368~370페이지에 설명되어 있습니다. 'eigen' 메서드는 테이퍼가 적용된 각 PSD 추정값에 이에 대응하는 슬레피안 테이퍼의 고유값(주파수 집중도)을 가중치로 적용합니다. 'unity' 메서드는 테이퍼가 적용된 각 PSD 추정값에 동일하게 가중치를 적용합니다.

e — DPSS(슬레피안) 시퀀스

행렬

DPSS(슬레피안) 시퀀스로, 행렬로 지정됩니다. x의 길이가 N이면 e는 N개 행을 가집니다. 행렬 e는 dpss의 출력값입니다.

v — DPSS(슬레피안) 시퀀스의 고유값

벡터

DPSS(슬레피안) 시퀀스의 고유값으로, 열 벡터로 지정됩니다. DPSS 시퀀스의 고유값은 시퀀스 에너지 비율이 분해능 대역폭 [–W, W]에 집중되어 있음을 나타냅니다. 고유값 범위는 구간 (0, 1) 내에 속하며, 대개 처음 2NW – 1개 고유값이 1에 가깝고 이후에는 0쪽으로 감소합니다. 벡터 v는 dpss의 출력값입니다.

dpss_params — dpss의 입력 인수

셀형 배열

dpss의 입력 인수로, 셀형 배열로 지정됩니다. dpss에 대한 첫 번째 입력 인수는 DPSS 시퀀스의 길이이며, x의 길이에서 구해지므로 dpss_params에서 생략됩니다.

예: pmtm(randn(1000,1),{2.5,3})은 시간과 반대역폭의 곱 2.5로 처음 3개의 슬레피안 시퀀스를 사용하여 랜덤 시퀀스의 PSD를 계산합니다.

freqrange — PSD 추정값의 주파수 범위

'onesided' | 'twosided' | 'centered'

PSD 추정값의 주파수 범위로, 'onesided', 'twosided' 또는 'centered'로 지정됩니다. 디폴트 값은 실수 값을 갖는 신호의 경우 'onesided'이고, 복소수 값을 갖는 신호의 경우 'twosided'입니다. 각 옵션에 대응되는 주파수 범위는 다음과 같습니다.

• 'onesided' — 실수 값을 갖는 입력 신호 x의 단측 PSD 추정값을 반환합니다. nfft가 짝수이면 pxx의 길이는 nfft/2 + 1이고 구간 [0,π] rad/sample에 대해 계산됩니다. nfft가 홀수이면 pxx의 길이는 (nfft + 1)/2이고 구간은 [0,π) rad/sample입니다. fs가 선택적으로 지정된 경우 이에 대응되는 구간은 짝수 길이와 홀수 길이 nfft 각각에 대해 [0,fs/2] cycles/unit time 및 [0,fs/2) cycles/unit time입니다.

  이 함수는 총 전력을 보존하기 위해 0 및 나이퀴스트 주파수를 제외한 모든 주파수에서 전력에 2를 곱합니다.

• 'twosided' — 실수 값 또는 복소수 값 입력 x의 양측 PSD 추정값을 반환합니다. 이 경우, pxx의 길이는 nfft이고 구간 [0,2π) rad/sample에 대해 계산됩니다. fs가 선택적으로 지정된 경우 구간은 [0,fs) cycles/unit time입니다.

• 'centered' — 실수 값 또는 복소수 값 입력 x에 대해 중심이 맞춰진 양측 PSD 추정값을 반환합니다. 이 경우, pxx의 길이는 nfft이며, nfft가 짝수 길이이면 구간 (–π,π] rad/sample에 대해 계산되고, nfft가 홀수 길이이면 구간 (–π,π) rad/sample에 대해 계산됩니다. fs가 선택적으로 지정된 경우 이에 대응되는 구간은 짝수 길이와 홀수 길이 nfft 각각에 대해 (–fs/2, fs/2] cycles/unit time 및 (–fs/2, fs/2) cycles/unit time입니다.

주파수로 구성된 벡터 w 또는 f와 freqrange를 입력값으로 지정하면 pmtm 함수는 freqrange를 무시합니다.

데이터형: char | string

probability — PSD 추정값의 신뢰구간

0.95 (디폴트 값) | (0,1) 범위의 스칼라

실제 PSD에 대한 포함 확률로, (0,1) 범위의 스칼라로 지정됩니다. 출력값 pxxc는 실제 PSD에 대한 probability × 100% 구간 추정값의 하한과 상한을 포함합니다.

출력 인수

pxx — PSD 추정값

벡터 | 행렬

PSD 추정값으로, 음이 아닌 실수 값 열 벡터나 행렬로 반환됩니다. pxx의 각 열은 x의 대응하는 열에 대한 PSD 추정값입니다. PSD 추정값의 단위는 단위 주파수당 시계열 데이터의 크기 단위의 제곱을 사용합니다. 예를 들어, 입력 데이터의 단위가 볼트이면 PSD 추정값의 단위는 단위 주파수당 제곱 볼트입니다. 볼트 단위의 시계열에 대해, 저항이 1Ω인 것으로 가정하고 샘플 레이트를 헤르츠로 지정할 경우 PSD 추정값의 단위는 헤르츠당 와트입니다.

w — 정규화 주파수

벡터

정규화 주파수로, 실수 값 열 벡터로 반환됩니다. pxx가 단측 PSD 추정값인 경우, w의 구간은 nfft가 짝수이면 [0,π]이고, nfft가 홀수이면 [0,π)입니다. pxx가 양측 PSD 추정값인 경우, w의 구간은 [0,2π)입니다. DC 기준 PSD 추정값의 경우, w의 구간은 nfft가 짝수이면 (–π,π]이고, nfft가 홀수이면 (–π,π)입니다.

f — 순환 주파수

벡터

순환 주파수로, 실수 값 열 벡터로 반환됩니다. 단측 PSD 추정값의 경우, f의 구간은 nfft가 짝수이면 [0,fs/2]이고, nfft가 홀수이면 [0,fs/2)입니다. 양측 PSD 추정값의 경우, f의 구간은 [0,fs)입니다. DC 기준 PSD 추정값의 경우, f의 구간은 nfft가 짝수 길이이면 (–fs/2, fs/2] cycles/unit time이고, nfft가 홀수 길이이면 (–fs/2, fs/2) cycles/unit time입니다.

pxxc — 신뢰한계

행렬

신뢰한계로, 실수 값 요소를 가진 행렬로 반환됩니다. 행렬의 행 크기는 PSD 추정값 pxx의 길이와 같습니다. pxxc는 pxx보다 두 배 더 많은 열을 가집니다. 홀수 열은 신뢰구간의 하한을 포함하고, 짝수 열은 상한을 포함합니다. 따라서 pxxc(m,2*n-1)은 추정값 pxx(m,n)에 대응하는 신뢰 하한이고 pxxc(m,2*n)은 추정값에 대응하는 신뢰 상한입니다. 신뢰구간의 포함 확률은 probability 입력값으로 결정됩니다.

데이터형: single | double

세부 정보

Thomson의 멀티테이퍼 스펙트럼 추정

주기도는 광의의 정상 과정의 실제 파워 스펙트럼 밀도(PSD)에 대한 일치 추정량(consistent estimator)이 아닙니다. 주기도의 가변성을 줄여서 일관된 PSD 추정값을 얻기 위해 멀티테이퍼 방법은 상호 직교 윈도우군, 즉 테이퍼로 얻은 수정된 주기도의 평균을 구합니다. 테이퍼는 상호 직교성 외에도 최적의 시간-주파수 집중도 속성을 가집니다. 테이퍼가 직교성과 시간-주파수 집중도를 모두 가진다는 사실은 멀티테이퍼 기법의 성패에 중요한 요인으로 작용합니다. Thomson의 멀티테이퍼 방법에 사용된 슬레피안 시퀀스에 대한 간략한 설명은 이산 장형 타원체(슬레피안) 시퀀스 항목을 참조하십시오.

멀티테이퍼 방법에서는 K개의 수정된 주기도를 사용합니다. 각각의 수정된 주기도는 서로 다른 슬레피안 시퀀스를 윈도우로 사용하여 얻은 것입니다.

$$S_k(f)=\Delta t{\left|{\displaystyle \sum _{n=0}^{N-1}{g}_k\left(n\right)}x\left(n\right)e^{-j2\pi fn\Delta t}\right|}^2$$

위 식은 k차 슬레피안 시퀀스 g_k(n)을 사용하여 얻은 수정된 주기도를 나타냅니다. 가장 단순한 형태의 멀티테이퍼 방법은 K개의 수정된 주기도에 대한 평균을 구하여 멀티테이퍼 PSD 추정값을 얻습니다.

$$S^{(\text{MT})}(f)=\frac{1}{K}{\displaystyle \sum _{k=0}^{K-1}{S}_k}(f).$$

[4]에서 소개된 Thomson의 멀티테이퍼 방법은 대략적으로 상관관계가 없는 PSD 추정값을 평균화한다는 점에서 Welch의 중첩 세그먼트 평균화 방법과 비슷합니다. 그러나 상관관계가 없는 PSD 추정값을 구하는 방식이 서로 다릅니다. 멀티테이퍼 방법은 수정된 주기도 각각에 담긴 전체 신호를 사용합니다. 슬레피안 테이퍼의 직교성을 통해, 다른 수정된 주기도와의 상관관계가 줄어듭니다. Welch의 방법은 수정된 주기도 각각의 신호 세그먼트를 사용하고, 분할을 통해 다른 수정된 주기도와의 상관을 줄입니다.

S^(MT)(f)에 대한 방정식은 pmtm의 'unity' 옵션에 해당합니다. 그러나 이산 장형 타원체(슬레피안) 시퀀스 항목에서 설명한 것처럼, 슬레피안 시퀀스는 원하는 주파수 대역에서 동일한 에너지 집중도를 갖지 않습니다. 슬레피안 시퀀스 차수가 높을수록, 집중도가 고유값으로 지정된 대역 [–W,W]에서의 시퀀스 에너지 집중도가 떨어집니다. 따라서 평균화하기 전에, K개의 수정된 주기도에 고유값을 사용하여 가중치를 적용하는 것이 좋습니다. 이는 pmtm에서의 'eigen' 옵션에 해당합니다.

수정된 주기도에 대해 가중 평균값을 구할 때 시퀀스 고유값을 사용하면 슬레피안 시퀀스의 주파수 집중도 속성이 고려됩니다. 그러나 확률 과정의 파워 스펙트럼 밀도와 슬레피안 시퀀스의 주파수 집중도 간의 상호 작용은 고려되지 않습니다. 특히 높은 차수의 슬레피안 시퀀스를 사용하는 경우, 수정된 주기도에서 전력이 매우 낮은 확률 과정의 주파수 영역이 안정적으로 추정되지 않습니다. 따라서 주파수 종속적 적응형 과정을 사용하는 것이 좋습니다. 이 과정은 슬레피안 시퀀스의 주파수 집중도뿐만 아니라 시계열의 전력 분포를 모두 고려합니다. 이 가변 가중치는 pmtm의 'adapt' 옵션에 해당하며, 멀티테이퍼 추정값을 계산하기 위한 디폴트 값입니다.

이산 장형 타원체(슬레피안) 시퀀스

슬레피안 시퀀스는 이산시간/연속 주파수 집중도 문제에서 유도됩니다. 인덱스가 0, 1, …, N – 1로 제한된 모든 ℓ² 시퀀스에 대해 이 문제는 |W| < 1/2Δt인 주파수 대역 [–W, W]에서 최대 에너지 집중도를 갖는 시퀀스를 찾습니다.

이는 N×N의 자기수반 양의 준정부호 연산자의 고유값과 이에 대응하는 고유벡터를 구하는 것과 같습니다. 따라서 고유값은 음이 아닌 실수이고, 중복되지 않은 고유값에 대응하는 고유벡터들은 서로 직교합니다. 이 특정 문제에서 고유값들은 1을 경계로 하며, 각 고유값은 주파수 구간 [–W, W]에서의 시퀀스 에너지 집중도 측정값에 해당합니다.

고유값 문제는 다음 방정식으로 지정됩니다.

$${\displaystyle \sum _{m=0}^{N-1}\frac{\sin \left(2\pi W\left(n-m\right)\right)}{\pi \left(n-m\right)}}{g}_k\left(m\right)={\lambda }_k\left(N,W\right)g_k\left(n\right),n,k=0,1,2,\dots ,N-1.$$

0차 DPSS 시퀀스 g₀은 가장 큰 고유값에 대응하는 고유벡터입니다. 1차 DPSS 시퀀스 g₁은 다음으로 가장 큰 고유값에 대응하는 고유벡터이고 0차 시퀀스와 직교합니다. 2차 DPSS 시퀀스 g₂는 세 번째로 큰 고유값에 대응하는 고유벡터이고 0차 및 1차 DPSS 시퀀스와 직교합니다. 연산자가 N×N이므로 N개 고유벡터가 있습니다. 그러나 지정된 시퀀스 길이 N과 지정된 대역폭 [–W, W]의 경우 고유값이 단위 값에 매우 가까운, 대략적으로 2NW – 1개의 DPSS 시퀀스가 있습니다. NW를 지정하려면 nw를 사용하십시오.

사인 테이퍼

슬레피안 시퀀스에 대한 대안으로 [3]에서 제안된 사인 테이퍼는 다음과 같이 정의됩니다.

$$g_k(n)=\sqrt{\frac{2}{N+1}}\sin \frac{\pi kn}{N+1},n,k=1,2,\dots ,N.$$

사인 테이퍼는 슬레피안 시퀀스와 달리 고유값 방정식을 설정하여 풀지 않고 직접 계산할 수 있습니다. 따라서 사인 테이퍼는 계산 속도가 훨씬 빠릅니다. 사인 테이퍼는 슬레피안 시퀀스와 비슷한 스펙트럼 집중도를 갖지만 스펙트럼 대역폭을 지정하기 위한 추가 파라미터가 필요하지 않습니다. 사인 테이퍼를 사용하여 계산된 PSD 추정값의 대역폭은 m을 사용하여 테이퍼의 개수를 변경함으로써 국소적으로 조정할 수 있습니다.

슬레피안 테이퍼와 사인 테이퍼 비교하기

시간과 반대역폭의 곱 3에 대응하는 처음 5개의 슬레피안 테이퍼를 생성합니다. 테이퍼 길이를 1000으로 지정합니다.

N = 1000;
nw = 3;
ns = 2*(nw)-1;

tprs = dpss(N,nw,ns);
lbs = "Slepian";

처음 5개의 사인 테이퍼를 생성합니다.

n = 1:N;
k = 1:ns;

tprs(:,:,2) = sqrt(2/(N+1))*sin(pi*n'*k/(N+1));
lbs(2) = "Sine";

두 개의 테이퍼 집합을 플로팅합니다.

for kj = 1:2
    subplot(2,1,kj)
    plot(tprs(:,:,kj))
    title(lbs(kj))
    legend(append('k = ',string(k+kj-2)), ...
        'Orientation','horizontal','Location','south')
    legend('boxoff')
    ylim([-0.09 0.07])
end

참고 문헌

[1] McCoy, Emma J., Andrew T. Walden, and Donald B. Percival. "Multitaper Spectral Estimation of Power Law Processes." IEEE^® Transactions on Signal Processing 46, no. 3 (March 1998): 655–68. https://doi.org/10.1109/78.661333.

[2] Percival, Donald B., and Andrew T. Walden. Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques. Cambridge; New York, NY, USA: Cambridge University Press, 1993.

[3] Riedel, Kurt S., and Alexander Sidorenko. “Minimum Bias Multiple Taper Spectral Estimation.” IEEE Transactions on Signal Processing 43, no. 1 (January 1995): 188–95. https://doi.org/10.1109/78.365298.

[4] Thomson, David J. "Spectrum estimation and harmonic analysis." Proceedings of the IEEE 70, no. 9 (1982): 1055–96. https://doi.org/10.1109/PROC.1982.12433.

확장 기능

C/C++ 코드 생성

MATLAB® Coder™를 사용하여 C 코드나 C++ 코드를 생성할 수 있습니다.

사용법 관련 참고 및 제한 사항:

• 입력 인수 nfft와 nw는 코드 생성 시 스칼라여야 합니다.

• nfft와 nw가 가변 크기인 경우 함수는 nfft를 주파수 벡터 f로 간주하고 nw를 슬레피안 테이퍼로 구성된 행렬 e로 간주합니다.

버전 내역

R2006a 이전에 개발됨

R2023b: C 및 C++ 코드 생성

pmtm 함수에 대한 C 및 C++ 코드를 생성합니다. 이 기능을 사용하려면 MATLAB^® Coder™가 있어야 합니다.

참고 항목

dpss | periodogram | pwelch

도움말 항목

• Bias and Variability in the Periodogram