Signal Processing Toolbox를 사용한 스펙트로그램 계산

Signal Processing Toolbox™는 비정상 신호의 스펙트로그램을 계산하는 3개의 함수를 제공합니다. 각 함수는 서로 다른 입력 인수, 디폴트 값, 출력값을 가집니다. 응용 분야에 따라 최적의 함수가 달라집니다.

비정상 신호는 시간에 따라 주파수 성분이 변하는 신호입니다. 단시간 푸리에 변환(STFT)은 신호가 변화함에 따라 이 주파수 성분이 어떻게 변하는지 분석하는 데 사용됩니다. STFT의 크기 제곱을 신호의 스펙트로그램 시간-주파수 표현이라고 합니다.

스펙트로그램은 여러 가능한 시간-주파수 표현 중 하나에 불과합니다. Signal Processing Toolbox와 Wavelet Toolbox™에서 사용 가능한 다른 시간-주파수 표현에 대한 개요는 Time-Frequency Gallery 항목을 참조하십시오. periodogram 함수를 사용하여 정상 신호를 처리하는 방법은 FFT를 사용한 파워 스펙트럼 밀도 추정값 항목을 참조하십시오.

스펙트로그램 계산 함수

Signal Processing Toolbox는 스펙트로그램 계산에 사용할 수 있는 다음과 같은 함수를 제공합니다.

spectrogram — 최대의 유연성을 위해 설계되었습니다. STFT와 세그먼트별 파워 스펙트럼 밀도 또는 파워 스펙트럼을 계산합니다. 재할당을 지원합니다.
stft — 가역성과 최대 제어를 위해 설계되었습니다. STFT를 계산합니다. dlstft와 stftLayer에서 사용됩니다. 다중채널 입력값을 지원합니다. 대응 함수인 istft 함수는 역 STFT를 계산합니다.
pspectrum — 사용 편의성을 위해 설계되었습니다. 파워 스펙트럼을 계산합니다. 신호 분석기에서 생성된 분석 스크립트에서 사용됩니다. 정상 신호 스펙트럼과 지속성 스펙트럼을 계산할 수 있습니다. 재할당을 지원합니다.

이후 설명에서는 spectrogram 함수를 기준으로 설명합니다.

범주	파라미터	함수
범주	파라미터	`spectrogram`	`stft`	`pspectrum`
입력값	신호	요소를 N_x개 가진 벡터	요소를 N_x개 가진 벡터 행을 N_x개 가진 행렬 행을 N_x개 가진 타임테이블	요소를 N_x개 가진 벡터 행을 N_x개 가진 타임테이블
	윈도우, g(n)	두 번째 위치 인수 (디폴트 값: 해밍 윈도우)	`Window` 이름-값 인수 (디폴트 값: 주기적 핸 윈도우)	카이저 윈도우만 해당
	윈도우 길이, M	샘플 수로 지정됨 (디폴트 값: ⌊N_x/4.5⌋)	샘플 수로 지정됨 (디폴트 값: 128)	`TimeResolution` 이름-값 인수
	누설, ℓ	윈도우에 따라 다름 카이저 윈도우를 사용하는 경우 β 형태 인자를 사용하여 조정	윈도우에 따라 다름 카이저 윈도우를 사용하는 경우 β 형태 인자를 사용하여 조정	카이저 윈도우 β 형태 인자와 관련된 `Leakage` 이름-값 인수: ℓ = 1 – β/40
	중첩, L	세 번째 위치 인수로 지정된 샘플 수 (디폴트 값: 윈도우 길이의 50%)	`OverlapLength` 이름-값 인수로 지정된 샘플 수 (디폴트 값: 윈도우 길이의 75%)	`OverlapPercent` 이름-값 인수로 지정된 세그먼트 길이의 백분율 (디폴트 값: $(1 - \frac{1}{2 \times ENBW - 1}) \times 100,$ , 여기서 ENBW는 윈도우의 등가 잡음 대역폭임)
	DFT 지점 개수, N_DFT	네 번째 위치 인수 (디폴트 값: max(256,2^{⌈log₂M⌉}))	`FFTLength` 이름-값 인수 (디폴트 값: 128)	항상 1024
	시간 정보	다섯 번째 위치 인수로 지정된 샘플 레이트	두 번째 위치 인수로 지정된 샘플 레이트 또는 시간 벡터	두 번째 위치 인수로 지정된 샘플 레이트 또는 시간 벡터

함수 호출	fs = 100; x = exp(2jpi20* ... (0:1/fs:2-1/fs)); M = 200; lk = 0.5; g = kaiser(M,40*(1-lk)); L = 100; Ndft = 1024;	sps = abs( ... spectrogram(x,g,L, ... Ndft,fs,"centered") ... ).^2;	sts = abs( ... stft(x,fs,Window=g, ... OverlapLength=L, ... FFTLength=Ndft) ... ).^2;	pss = pspectrum(x,fs, ... "spectrogram", ... TimeResolution=M/fs, ... Leakage=lk, ... OverlapPercent=L/M100 ... )sum(g)^2;

편의 플롯	fs = 6e2; ts = 0:1/fs:2.05; x = vco(sin(2pits).* ... exp(-ts),[0.1 0.4]fs,fs); M = 32; lk = 0.9; g = kaiser(M,40(1-lk)); L = 22; Ndft = 1024; `spectrogram`의 경우 `"power","yaxis"` 추가 `stft`의 경우 `FrequencyRange="onesided"` 추가

출력값	주파수 범위	`freqrange` 인수를 사용하여 제어됨: `"onesided"` — N_DFT 값이 짝수인 경우 주파수 구간은 영주파수와 나이퀴스트 주파수 f_s/2에서 닫힙니다. N_DFT 값이 홀수인 경우 주파수 구간은 영주파수에서 닫히고 f_s/2에서 열립니다. (실수 값 신호의 디폴트 값.) `"twosided"` — 모든 N_DFT 값에 대해 주파수 구간은 영주파수에서 닫히고 f_s에서 열립니다. `"centered"`— N_DFT 값이 짝수인 경우 주파수 구간은 –f_s/2에서 열리고 f_s/2에서 닫힙니다. N_DFT 값이 홀수인 경우 주파수 구간은 양끝에서 열립니다. (복소수 값 신호의 디폴트 값.) 사용자는 STFT와 스펙트로그램을 계산할 주파수 벡터를 지정할 수 있습니다.	`FrequencyRange` 이름-값 인수를 사용하여 제어됩니다. `"onesided"` — `spectrogram`과 동일합니다. `"twosided"` — `spectrogram`과 동일합니다. `"centered"` — `spectrogram`과 동일합니다. (실수 값 신호 및 복소수 값 신호의 디폴트 값.)	`TwoSided` 논리형 이름-값 인수를 사용하여 제어됩니다. `false` — 영주파수와 f_s/2에서 구간이 닫힙니다. (실수 값 신호의 디폴트 값.) `true` – –f_s/2와 f_s/2에서 구간이 닫힙니다. (복소수 값 신호의 디폴트 값.)
	시간 간격	전체 윈도우가 적용된 마지막 세그먼트 이후의 신호는 잘립니다. 시간 값은 세그먼트의 중심 지점입니다.	전체 윈도우가 적용된 마지막 세그먼트 이후의 신호는 잘립니다. 시간 값은 세그먼트의 중심 지점입니다.	전체 윈도우가 적용된 마지막 세그먼트 이후의 신호는 0으로 채워집니다. 시간 값은 세그먼트의 중심 지점입니다.
	정규화	첫 번째 출력 인수는 STFT입니다. 그 크기를 제곱한 것이 스펙트로그램입니다. 네 번째 출력 인수는 크기 제곱입니다. 스펙트로그램을 구하려면 (Σ_ng(n))²을 곱하고 `"power"` 옵션을 지정합니다.	첫 번째 출력 인수는 STFT입니다. 그 크기를 제곱한 것이 스펙트로그램입니다.	첫 번째 출력 인수는 크기 제곱입니다. 스펙트로그램을 구하려면 (Σ_ng(n))²을 곱합니다.
	PSD 및 파워 스펙트럼	`spectrogram`의 네 번째 출력 인수에는 세그먼트 파워 스펙트럼 밀도 또는 세그먼트 파워 스펙트럼이 포함됩니다. 스펙트로그램은 파워 스펙트럼에 윈도우 요소 합계의 제곱을 곱한 값입니다. `spectrumtype` 인수: `"psd"`— ENBW를 곱하여 파워 스펙트럼을 구함 (디폴트 값) `"power"` — 파워 스펙트럼	출력 인수는 STFT입니다.	출력 인수는 파워 스펙트럼입니다. PSD를 구하려면 ENBW로 나눕니다.

예제		스펙트로그램의 디폴트 값 spectrogram 함수와 STFT 정의 비교하기 오디오 신호의 처프(Chirp) 추적하기	단시간 푸리에 변환 spectrogram 함수와 stft 함수 비교하기 다중채널 신호의 STFT	정현파의 파워 스펙트럼 spectrogram 함수와 pspectrum 함수 비교하기 음악 신호에 대역저지 필터 적용하기

STFT 및 스펙트로그램 정의

신호의 STFT는 길이가 M인 분석 윈도우 g(n)를 신호에 대해 적용한 후, 윈도우가 적용된 데이터의 이산 푸리에 변환을 계산하는 방식으로 계산됩니다. 이 윈도우는 원래 신호에 대해 R개 샘플의 간격만큼 건너뛰며 적용되며, 이는 인접 세그먼트 간의 샘플 중첩 L = M – R과 동일합니다. 대부분 윈도우 함수는 스펙트럼 링잉 현상을 피하기 위해 가장자리에서 값이 감소합니다. 윈도우가 적용된 각 세그먼트의 DFT는 각 시간 및 주파수 점에 대한 크기와 위상이 포함된 복소수 값 행렬에 추가됩니다. STFT 행렬에는 다음 개수의 열이 있습니다.

$k = ⌊ \frac{N_{x} - L}{M - L} ⌋$

여기서 N_x는 신호 x(n)의 길이이고 ⌊⌋ 기호는 바닥 함수를 나타냅니다. 행렬의 행 개수는 중심이 맞춰진 양측 변환의 경우 DFT 지점 개수 N_DFT와 같고, 실수 값 신호의 단측 변환의 경우 N_DFT/2에 가까운 홀수입니다.

STFT 행렬 $X (f) = [\begin{matrix} X_{1} (f) & X_{2} (f) & X_{3} (f) & \dots & X_{k} (f) \end{matrix}]$ 의 m번째 열에는 시간 mR을 중심으로 하는 윈도우가 적용된 데이터의 DFT가 포함되어 있습니다.

$X_{m} (f) = \sum_{n = - \infty}^{\infty} x (n) g (n - m R) e^{- j 2 π f n} .$

`spectrogram` 함수와 STFT 정의 비교하기

라이브 스크립트 열기

2초 동안 600Hz로 샘플링된 복소수 값의 볼록 2차 처프로 구성된 신호를 생성합니다. 처프의 초기 주파수는 250Hz이고 최종 주파수는 50Hz입니다.

fs = 6e2;
ts = 0:1/fs:2;
x = chirp(ts,250,ts(end),50,"quadratic",0,"convex","complex");

spectrogram 함수

spectrogram 함수를 사용하여 신호의 STFT를 계산합니다.

신호를 각 샘플 길이가 $M = 49$ 인 세그먼트로 나눕니다.
인접 세그먼트 간에 중첩되는 샘플을 $L = 11$ 로 지정합니다.
더 짧은 최종 세그먼트는 무시합니다.
바틀렛 윈도우로 각 세그먼트에 윈도우를 적용합니다.
$N_{DFT} = 1024$ 개 점에서 각 세그먼트의 이산 푸리에 변환을 평가합니다. 기본적으로 spectrogram은 복소수 값 신호에 대해 양측 변환을 계산합니다.

M = 49;
L = 11;
g = bartlett(M);
Ndft = 1024;

[s,f,t] = spectrogram(x,g,L,Ndft,fs);

waterplot 함수를 사용하여 STFT의 크기 제곱으로 정의되는 스펙트로그램을 계산하고 표시합니다.

waterplot(s,f,t)

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

STFT 정의

정의를 사용하여 $N_{x}$ 개 샘플 신호의 STFT를 계산합니다. 신호를 $⌊ \frac{N_{x} - L}{M - L} ⌋$ 중첩 세그먼트로 나눕니다. 각 세그먼트에 윈도우를 적용하고 $N_{DFT}$ 개 점에서 이산 푸리에 변환을 평가합니다.

segs = framesig(1:length(x),M,OverlapLength=L);
X = fft(x(segs).*g,Ndft);

STFT의 시간과 주파수 범위를 계산합니다.

시간 값을 찾으려면 시간 벡터를 중첩 세그먼트로 나눕니다. 시간 값은 세그먼트의 중간점이며 각 세그먼트는 하한에서 열린 구간으로 처리됩니다.
주파수 값을 찾으려면 영주파수에서 닫히고 하한에서 열리는 나이퀴스트 구간을 지정합니다.

framedT = ts(segs);
tint = mean(framedT(2:end,:));

fint = 0:fs/Ndft:fs-fs/Ndft;

spectrogram의 출력값을 정의와 비교합니다. 스펙트로그램을 표시합니다.

maxdiff = max(max(abs(s-X)))

maxdiff = 
0

waterplot(X,fint,tint)

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency (Hz)")
    ylabel("Time (s)")
end

`spectrogram` 함수와 `stft` 함수 비교하기

라이브 스크립트 열기

2초 동안 1.4kHz로 샘플링된 처프로 구성된 신호를 생성합니다. 측정 시간 동안 처프의 주파수는 600Hz에서 100Hz로 선형적으로 감소합니다.

fs = 1400;
x = chirp(0:1/fs:2,600,2,100);

stft 디폴트 값

spectrogram과 stft 함수를 사용하여 신호의 STFT를 계산합니다. stft 함수의 디폴트 값을 사용합니다.

신호를 한 세그먼트당 128개 샘플이 있도록 여러 세그먼트로 나누고 주기적 핸 윈도우를 사용하여 각 세그먼트에 윈도우를 적용합니다.
인접 세그먼트 간에 96개 샘플이 중첩되도록 지정합니다. 이 길이는 윈도우 길이의 75%에 해당합니다.
128개의 DFT 점을 지정하고 헤르츠 단위로 표현되는 STFT를 영주파수의 중심에 둡니다.

두 결과가 동일한지 확인합니다.

M = 128;
g = hann(M,"periodic");
L = 0.75*M;
Ndft = 128;

[sp,fp,tp] = spectrogram(x,g,L,Ndft,fs,"centered");

[s,f,t] = stft(x,fs);

dff = max(max(abs(sp-s)))

dff = 
0

mesh 함수를 사용하여 두 출력값을 플로팅합니다.

nexttile
mesh(tp,fp,abs(sp).^2)
title("spectrogram")
view(2), axis tight

nexttile
mesh(t,f,abs(s).^2)
title("stft")
view(2), axis tight

Figure contains 2 axes objects. Axes object 1 with title spectrogram contains an object of type surface. Axes object 2 with title stft contains an object of type surface.

spectrogram 디폴트 값

spectrogram 함수의 디폴트 값을 사용하여 계산을 반복합니다.

신호를 길이가 $M = ⌊ N_{x} / 4.5 ⌋$ 인 세그먼트로 나눕니다. 여기서 $N_{x}$ 는 신호의 길이입니다. 해밍 윈도우로 각 세그먼트에 윈도우를 적용합니다.
세그먼트 간에 50% 중첩을 지정합니다.
$\max (256, 2^{⌈ \log_{2} M ⌉})$ 개의 점을 사용하여 FFT를 계산합니다. 양의 정규화 주파수에 대해서만 스펙트로그램을 계산합니다.

M = floor(length(x)/4.5);
g = hamming(M);
L = floor(M/2);
Ndft = max(256,2^nextpow2(M));

[sx,fx,tx] = spectrogram(x);

[st,ft,tt] = stft(x,Window=g,OverlapLength=L, ...
    FFTLength=Ndft,FrequencyRange="onesided");

dff = max(max(sx-st))

dff = 
0

waterplot 함수를 사용하여 두 출력값을 플로팅합니다. 두 경우 모두 주파수 축을 $π$ 로 나눕니다. stft 출력값의 경우 샘플 개수를 유효 샘플 레이트 $2 π$ 로 나눕니다.

figure
nexttile
waterplot(sx,fx/pi,tx)
title("spectrogram")

nexttile
waterplot(st,ft/pi,tt/(2*pi))
title("stft")

$Figure contains 2 axes objects. Axes object 1 with title spectrogram, xlabel Frequency/\pi, ylabel Samples contains an object of type patch. Axes object 2 with title stft, xlabel Frequency/\pi, ylabel Samples contains an object of type patch.$

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency/\pi")
    ylabel("Samples")
end

`spectrogram` 함수와 `pspectrum` 함수 비교하기

라이브 스크립트 열기

전압 제어 발진기와 3개의 가우스 원자로 구성된 신호를 생성합니다. 신호는 2초 동안 $f_{s} = 2$ kHz로 샘플링됩니다.

fs = 2000;
tx = 0:1/fs:2;
gaussFun = @(A,x,mu,f) exp(-(x-mu).^2/(2*0.03^2)).*sin(2*pi*f.*x)*A';
s = gaussFun([1 1 1],tx',[0.1 0.65 1],[2 6 2]*100)*1.5;
x = vco(chirp(tx+.1,0,tx(end),3).*exp(-2*(tx-1).^2),[0.1 0.4]*fs,fs);
x = s+x';

단시간 푸리에 변환

pspectrum 함수를 사용하여 STFT를 계산합니다.

$N_{x}$ 개 샘플 신호를 $80 / 2000 = 40$ 밀리초의 시간 분해능에 상응하는 샘플 길이 $M = 80$ 인 세그먼트로 나눕니다.
인접 세그먼트 간 중첩 샘플 개수 $L = 16$ 을 지정하거나 20%가 중첩되도록 지정합니다.
카이저 윈도우로 각 세그먼트에 윈도우를 적용하고 누설 $ℓ = 0.7$ 을 지정합니다.

M = 80;
L = 16;
lk = 0.7;

[S,F,T] = pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk);

spectrogram 함수로 구한 결과와 비교합니다.

윈도우 길이를 지정하고 샘플에서 직접 중첩합니다.
pspectrum은 항상 카이저 윈도우를 $g (n)$ 으로 사용합니다. 누설 $ℓ$ 과 윈도우의 형태 인자 $β$ 는 $β = 40 \times (1 - ℓ)$ 의 관계를 갖습니다.
pspectrum은 이산 푸리에 변환을 계산할 때 항상 $N_{DFT} = 1024$ 개 점을 사용합니다. 양측 주파수 범위 또는 중심 주파수 범위에서 변환을 계산하려는 경우 이 숫자를 지정할 수 있습니다. 그러나 실수 신호의 디폴트인 단측 변환의 경우 spectrogram은 $1024 / 2 + 1 = 513$ 개 점을 사용합니다. 또는 이 예제와 같이 변환을 계산할 주파수 벡터를 지정할 수 있습니다.
신호를 정확히 $k = ⌊ \frac{N_{x} - L}{M - L} ⌋$ 개 세그먼트로 나눌 수 없는 경우 spectrogram은 신호를 자르고 pspectrum은 신호에 0을 채워 여분의 세그먼트를 만듭니다. 출력값을 동등하게 만들려면 마지막 세그먼트와 시간 벡터의 마지막 요소를 제거하십시오.
spectrogram은 크기 제곱이 스펙트로그램인 STFT를 반환합니다. pspectrum은 이미 제곱되었지만 제곱 전에 인자 $\sum_{n} g (n)$ 으로 나눈 세그먼트별 파워 스펙트럼을 반환합니다.
단측 변환의 경우 pspectrum은 스펙트로그램에 추가 인자 2를 추가합니다.

g = kaiser(M,40*(1-lk));

k = (length(x)-L)/(M-L);
if k~=floor(k)
    S = S(:,1:floor(k));
    T = T(1:floor(k));
end

[s,f,t] = spectrogram(x/sum(g)*sqrt(2),g,L,F,fs);

waterplot 함수를 사용하여 두 함수로 계산된 스펙트로그램을 표시합니다.

subplot(2,1,1)
waterplot(sqrt(S),F,T)
title("pspectrum")

subplot(2,1,2)
waterplot(s,f,t)
title("spectrogram")

Figure contains 2 axes objects. Axes object 1 with title pspectrum, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title spectrogram, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

maxd = max(max(abs(abs(s).^2-S)))

maxd = 
2.4419e-08

파워 스펙트럼 및 편의 플롯

spectrogram 함수에는 세그먼트별 파워 스펙트럼 또는 파워 스펙트럼 밀도에 대응하는 네 번째 인수가 있습니다. ps 인수는 pspectrum의 출력값과 유사하게 이미 제곱되었으며 정규화 인자 $\sum_{n} g (n)$ 을 포함합니다. 실수 신호의 단측 스펙트로그램의 경우 추가 인자 2를 꼭 포함해야 합니다. 함수의 스케일링 인수를 "power"로 설정합니다.

[~,~,~,ps] = spectrogram(x*sqrt(2),g,L,F,fs,"power");

max(abs(S(:)-ps(:)))

ans = 
2.4419e-08

출력 인수 없이 호출하면 pspectrum과 spectrogram 모두 데시벨 단위로 신호의 스펙트로그램을 플로팅합니다. 단측 스펙트로그램의 경우 인자 2를 포함시킵니다. 두 플롯의 컬러맵을 동일하게 설정합니다. x 제한을 동일한 값으로 설정하여 pspectrum 플롯 끝에 여분의 세그먼트가 보이도록 합니다. spectrogram 플롯에서 y축에 주파수를 표시합니다.

subplot(2,1,1)
pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk)
title("pspectrum")
cc = clim;
xl = xlim;

subplot(2,1,2)
spectrogram(x*sqrt(2),g,L,F,fs,"power","yaxis")
title("spectrogram")
clim(cc)
xlim(xl)

Figure contains 2 axes objects. Axes object 1 with title pspectrum, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image. Axes object 2 with title spectrogram, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image.

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency (Hz)")
    ylabel("Time (s)")
end

중심 스펙트로그램과 단측 스펙트로그램 계산하기

라이브 스크립트 열기

2초 동안 2kHz로 샘플링된 실수 값 처프로 구성된 신호를 생성합니다.

fs = 2000;
tx = 0:1/fs:2;
x = vco(-chirp(tx,0,tx(end),2).*exp(-3*(tx-1).^2), ...
    [0.1 0.4]*fs,fs).*hann(length(tx))';

양측 스펙트로그램

신호의 양측 STFT를 계산하고 플로팅합니다.

신호를 각 샘플 길이가 $M = 73$ 인 세그먼트로 나눕니다.
인접 세그먼트 간에 중첩되는 샘플을 $L = 24$ 로 지정합니다.
더 짧은 최종 세그먼트는 무시합니다.
플랫 탑(Flat Top) 윈도우로 각 세그먼트에 윈도우를 적용합니다.
$N_{DFT} = 895$ 개 점에서 각 세그먼트의 이산 푸리에 변환을 평가합니다. 이때 홀수 개수라는 점에 유의하십시오.

M = 73;
L = 24;
g = flattopwin(M);
Ndft = 895;
neven = ~mod(Ndft,2);

[stwo,f,t] = spectrogram(x,g,L,Ndft,fs,"twosided");

출력 인수 없이 spectrogram 함수를 사용하여 양측 스펙트로그램을 플로팅합니다.

spectrogram(x,g,L,Ndft,fs,"twosided","power","yaxis")

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Frequency (kHz) contains an object of type image.

정의를 사용하여 양측 스펙트로그램을 계산합니다. 세그먼트의 샘플 길이가 $M$ 이고 인접 세그먼트 간의 중첩이 $L$ 이 되도록 신호를 나눕니다. 각 세그먼트에 윈도우를 적용하고 $N_{DFT}$ 개 점에서 이산 푸리에 변환을 계산합니다.

y = framesig(x,M,Window=g,OverlapLength=L);
Xtwo = fft(y,Ndft);

시간과 주파수 범위를 계산합니다.

시간 값을 찾으려면 시간 벡터를 중첩 세그먼트로 나눕니다. 시간 값은 세그먼트의 중간점이며 각 세그먼트는 하한에서 열린 구간으로 처리됩니다.
주파수 값을 찾으려면 영주파수에서 닫히고 상한에서 열리는 나이퀴스트 구간을 지정합니다.

framedT = framesig(tx,M,OverlapLength=L);
ttwo = mean(framedT(2:end,:));

ftwo = 0:fs/Ndft:fs*(1-1/Ndft);

spectrogram의 출력값을 정의와 비교합니다. waterplot 함수를 사용하여 스펙트로그램을 표시합니다.

diffs = [max(max(abs(stwo-Xtwo)));
    max(abs(f-ftwo'));
    max(abs(t-ttwo))]

diffs = 3×1
10^-12 ×

         0
    0.2274
    0.0002

figure
nexttile
waterplot(Xtwo,ftwo,ttwo)
title("Two-Sided, Definition")

nexttile
waterplot(stwo,f,t)
title("Two-Sided, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title Two-Sided, Definition, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title Two-Sided, spectrogram Function, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

중심 스펙트로그램

신호의 중심 스펙트로그램을 계산합니다.

양측 STFT에 사용한 것과 같은 시간 값을 사용합니다.
fftshift 함수를 사용하여 STFT의 영주파수 성분을 스펙트럼의 중심으로 이동합니다.
홀수 값 $N_{DFT}$ 의 경우 주파수 구간이 양끝에서 열립니다. 짝수 값 $N_{DFT}$ 의 경우 주파수 구간이 하한에서 열리고 상한에서 닫힙니다.

출력값을 비교하고 스펙트로그램을 표시합니다.

tcen = ttwo;

if ~neven
    Xcen = fftshift(Xtwo,1);
    fcen = -fs/2*(1-1/Ndft):fs/Ndft:fs/2;
else
    Xcen = fftshift(circshift(Xtwo,-1),1);
    fcen = (-fs/2*(1-1/Ndft):fs/Ndft:fs/2)+fs/Ndft/2;
end

[scen,f,t] = spectrogram(x,g,L,Ndft,fs,"centered");

diffs = [max(max(abs(scen-Xcen)));
    max(abs(f-fcen'));
    max(abs(t-tcen))]

diffs = 3×1
10^-12 ×

         0
    0.2274
    0.0002

figure
nexttile
waterplot(Xcen,fcen,tcen)
title("Centered, Definition")

nexttile
waterplot(scen,f,t)
title("Centered, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title Centered, Definition, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title Centered, spectrogram Function, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

단측 스펙트로그램

신호의 단측 스펙트로그램을 계산합니다.

양측 STFT에 사용한 것과 같은 시간 값을 사용합니다.
홀수 값 $N_{DFT}$ 의 경우 단측 STFT는 양측 STFT의 처음 $(N_{DFT} + 1) / 2$ 개 행으로 구성됩니다. 짝수 값 $N_{DFT}$ 의 경우 단측 STFT는 양측 STFT의 처음 $N_{DFT} / 2 + 1$ 개 행으로 구성됩니다.
홀수 값 $N_{DFT}$ 의 경우 주파수 구간은 영주파수에서 닫히고 나이퀴스트 주파수에서 열립니다. 짝수 값 $N_{DFT}$ 의 경우 주파수 구간은 양끝에서 닫힙니다.

출력값을 비교하고 스펙트로그램을 표시합니다. 실수 값 신호의 경우 "onesided" 인수는 선택 사항입니다.

tone = ttwo;

if ~neven
    Xone = Xtwo(1:(Ndft+1)/2,:);
else
    Xone = Xtwo(1:Ndft/2+1,:);
end

fone = 0:fs/Ndft:fs/2;

[sone,f,t] = spectrogram(x,g,L,Ndft,fs);

diffs = [max(max(abs(sone-Xone)));
    max(abs(f-fone'));
    max(abs(t-tone))]

diffs = 3×1
10^-12 ×

         0
    0.1137
    0.0002

figure
nexttile
waterplot(Xone,fone,tone)
title("One-Sided, Definition")

nexttile
waterplot(sone,f,t)
title("One-Sided, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title One-Sided, Definition, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title One-Sided, spectrogram Function, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

function waterplot(s,f,t)
% Waterfall plot of spectrogram
waterfall(f,t,abs(s)'.^2)
set(gca,XDir="reverse",View=[30 50])
xlabel("Frequency (Hz)")
ylabel("Time (s)")
end

세그먼트 PSD와 파워 스펙트럼 계산하기

라이브 스크립트 열기

spectrogram 함수는 각 세그먼트의 파워 스펙트럼 밀도(PSD) 또는 파워 스펙트럼의 행렬을 네 번째 출력 인수로 가집니다. 파워 스펙트럼은 PSD에 윈도우의 등가 잡음 대역폭(ENBW)을 곱한 것과 같습니다.

1초 동안 1kHz로 샘플링된 로그 처프로 구성된 신호를 생성합니다. 처프는 초기 주파수가 400Hz이며 측정이 끝나는 시점에는 10Hz까지 감소합니다.

fs = 1000;
tt = 0:1/fs:1-1/fs;
y = chirp(tt,400,tt(end),10,"logarithmic");

샘플 레이트를 사용하여 PSD와 파워 스펙트럼 분할하기

신호를 한 세그먼트당 102개 샘플이 있도록 여러 세그먼트로 나누고 핸 윈도우를 사용하여 각 세그먼트에 윈도우를 적용합니다. 인접 세그먼트 간에 12개 샘플이 중첩되도록 지정하고 1024개 DFT 점을 지정합니다.

M = 102;
g = hann(M);
L = 12;
Ndft = 1024;

디폴트 PSD 스펙트럼 유형을 사용하여 신호의 스펙트로그램을 계산합니다. STFT와 세그먼트 파워 스펙트럼 밀도의 배열을 출력합니다.

[s,f,t,p] = spectrogram(y,g,L,Ndft,fs);

스펙트럼 유형을 "power"로 지정하여 계산을 반복합니다. STFT와 세그먼트 파워 스펙트럼의 배열을 출력합니다.

[r,~,~,q] = spectrogram(y,g,L,Ndft,fs,"power");

두 경우 모두 스펙트로그램이 동일한지 확인합니다. 주파수에 로그 스케일을 사용하여 스펙트로그램을 플로팅합니다.

max(max(abs(s).^2-abs(r).^2))

ans = 
0

waterfall(f,t,abs(s)'.^2)
set(gca,XScale="log",...
    XDir="reverse",View=[30 50])

Figure contains an axes object. The axes object contains an object of type patch.

파워 스펙트럼이 파워 스펙트럼 밀도에 윈도우의 ENBW를 곱한 것과 같은지 확인합니다.

max(max(abs(q-p*enbw(g,fs))))

ans = 
1.1102e-16

세그먼트 파워 스펙트럼의 행렬이 스펙트로그램에 비례하는지 확인합니다. 비례 인자는 윈도우 요소 합계의 제곱입니다.

max(max(abs(s).^2-q*sum(g)^2))

ans = 
0

정규화 주파수를 사용하여 PSD와 파워 스펙트럼 분할하기

계산을 반복하되 이번에는 정규화 주파수를 사용합니다. 샘플 레이트를 $2 π$ 로 지정하면 결과가 같아집니다.

[~,~,~,pn] = spectrogram(y,g,L,Ndft);
[~,~,~,qn] = spectrogram(y,g,L,Ndft,"power");

max(max(abs(qn-pn*enbw(g,2*pi))))

ans = 
1.1102e-16

참고 항목

앱

신호 분석기

함수

dlstft | dlistft | pspectrum | spectrogram | stft | istft | xspectrogram