stft

단시간 푸리에 변환

구문

s = stft(x)

s = stft(x,fs)

s = stft(x,ts)

s = stft(___,Name=Value)

[s,f] = stft(___)

[s,f,t] = stft(___)

stft(___)

설명

s = stft(x)는 x의단시간 푸리에 변환(STFT)을 반환합니다.

s = stft(x,fs)는 샘플 레이트 fs를 사용하여 x의 STFT를 반환합니다.

s = stft(x,ts)는 샘플 시간 ts를 사용하여 x의 STFT를 반환합니다.

s = stft(___,Name=Value)는 이름-값 인수를 사용하여 추가 옵션을 지정합니다. 이런 옵션은 FFT 윈도우와 길이가 있습니다. 이 인수는 위에 열거된 입력 구문에 추가할 수 있습니다.

[s,f] = stft(___)는 STFT가 계산되는 주파수 f를 반환합니다.

[s,f,t] = stft(___)는 STFT가 계산되는 시간을 반환합니다.

stft(___)에 출력 인수를 지정하지 않으면 현재 Figure 창에 STFT의 크기 제곱(단위: 데시벨)을 플로팅합니다.

예제

단시간 푸리에 변환

정현적으로 변하는 주파수를 갖는 처프를 생성합니다. 신호는 2초 동안 10kHz로 샘플링됩니다.

fs = 10e3;
t = 0:1/fs:2;
x = vco(sin(2*pi*t),[0.1 0.4]*fs,fs);

처프에 대해 단시간 푸리에 변환을 계산합니다. 신호를 한 세그먼트당 256개 샘플이 있도록 여러 세그먼트로 나누고 형태 파라미터 β = 5인 카이저 윈도우를 사용하여 각 세그먼트에 윈도우를 적용합니다. 인접 세그먼트 간에 220개 샘플이 중첩되도록 지정하고 DFT 길이를 512로 지정합니다. STFT가 계산되는 주파수와 시간 값을 출력합니다.

[s,f,t] = stft(x,fs,Window=kaiser(256,5),OverlapLength=220,FFTLength=512);

STFT의 크기 제곱을 스펙트로그램이라고도 합니다. 스펙트로그램(단위: 데시벨)을 플로팅합니다. 범위가 60dB이고 마지막 색이 스펙트로그램의 최댓값에 대응하도록 컬러맵을 지정합니다.

sdb = mag2db(abs(s));
mesh(t,f/1000,sdb);

cc = max(sdb(:))+[-60 0];
ax = gca;
ax.CLim = cc;
view(2)
colorbar

출력 인수 없이 stft 함수를 호출하여 동일한 플롯을 얻습니다.

stft(x,fs,Window=kaiser(256,5),OverlapLength=220,FFTLength=512)

2차 처프의 STFT

2초 동안 1kHz로 샘플링된 2차 처프(Chirp)를 생성합니다. 순시 주파수는 $\mathit{t}=0$초에서 100Hz, $\mathit{t}=1$초에서 200Hz를 넘어섭니다.

ts = 0:1/1e3:2;

f0 = 100;
f1 = 200;

x = chirp(ts,f0,1,f1,"quadratic",[],"concave");

지속 시간이 1ms인 2차 처프의 STFT를 계산하고 표시합니다.

d = seconds(1e-3);
win = hamming(100,"periodic");

stft(x,d,Window=win,OverlapLength=98,FFTLength=128);

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

spectrogram 디폴트 값

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

• 신호를 길이가 $\mathit{M}=\lfloor {\mathit{N}}_{\mathit{x}}/4.5\rfloor$인 세그먼트로 나눕니다. 여기서 ${\mathit{N}}_{\mathit{x}}$는 신호의 길이입니다. 해밍 윈도우로 각 세그먼트에 윈도우를 적용합니다.

• 세그먼트 간에 50% 중첩을 지정합니다.

• $\max \left(256,2^{\lceil {\log }_2\mathit{M}\rceil }\right)$개의 점을 사용하여 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 함수를 사용하여 두 출력값을 플로팅합니다. 두 경우 모두 주파수 축을 $\pi$로 나눕니다. stft 출력값의 경우 샘플 개수를 유효 샘플 레이트 $2\pi$로 나눕니다.

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

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

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

STFT 주파수 범위

4초 동안 5kHz로 샘플링된 신호를 생성합니다. 신호는 증가하는 추세를 갖는 변동적인 주파수와 진동하는 진폭이 있는 영역으로 구분된, 지속 기간이 줄어드는 펄스 세트로 구성됩니다. 신호를 플로팅합니다.

fs = 5000;
t = 0:1/fs:4-1/fs;

x = besselj(0,600*(sin(2*pi*(t+1).^3/30).^5));

plot(t,x)

신호의 단측, 양측 및 중심이 맞춰진 단시간 푸리에 변환을 계산합니다. 모든 경우에, 형태 인자 $\beta =10$을 갖는 202개 샘플로 구성된 카이저 윈도우를 사용하여 신호 세그먼트에 윈도우를 적용합니다. 각 변환을 계산하는 데 사용된 주파수 범위를 표시합니다.

rngs = ["onesided" "twosided" "centered"];

for kj = 1:length(rngs)
    
    opts = {'Window',kaiser(202,10),'FrequencyRange',rngs(kj)};

    [~,f] = stft(x,fs,opts{:});
    subplot(length(rngs),1,kj)
    stft(x,fs,opts{:})
    title(sprintf('''%s'': [%5.3f, %5.3f] kHz',rngs(kj),[f(1) f(end)]/1000))

end

계산을 반복하되, 이번에는 카이저 윈도우의 길이를 홀수인 203으로 변경합니다. 'twosided' 주파수 구간은 변경되지 않습니다. 나머지 두 개의 주파수 구간은 높은 쪽 끝에서 열린 상태가 됩니다.

for kj = 1:length(rngs)
    
    opts = {'Window',kaiser(203,10),'FrequencyRange',rngs(kj)};

    [~,f] = stft(x,fs,opts{:});
    subplot(length(rngs),1,kj)
    stft(x,fs,opts{:})
    title(sprintf('''%s'': [%5.3f, %5.3f] kHz',rngs(kj),[f(1) f(end)]/1000))

end

다중채널 신호의 STFT

1초 동안 1kHz로 샘플링된 3개의 서로 다른 처프로 구성된 3채널 신호를 생성합니다.

1. 첫 번째 채널은 오목 2차 처프로 구성되며, t = 0 지점에서는 순시 주파수가 100Hz이고 t = 1 지점에서는 300Hz를 넘어섭니다. 초기 위상이 45도입니다.

2. 두 번째 채널은 볼록 2차 처프로 구성되며, t = 0 지점에서는 순시 주파수가 100Hz이고 t = 1 지점에서는 500Hz를 넘어섭니다.

3. 두 번째 채널은 로그 2차 처프로 구성되며, t = 0 지점에서는 순시 주파수가 300Hz이고 t = 1 지점에서는 500Hz를 넘어섭니다.

길이가 128인 주기 해밍 윈도우 및 50개 샘플의 중첩 길이를 사용하여 다중채널 신호의 STFT를 계산합니다.

fs = 1e3;
t = 0:1/fs:1-1/fs;

x = [chirp(t,100,1,300,'quadratic',45,'concave');
      chirp(t,100,1,500,'quadratic',[],'convex');
      chirp(t,300,1,500,'logarithmic')]'; 
  
[S,F,T] = stft(x,fs,'Window',hamming(128,'periodic'),'OverlapLength',50);

각 채널의 STFT를 폭포 플롯으로 시각화합니다. 헬퍼 함수 helperGraphicsOpt를 사용하여 좌표축의 동작을 제어합니다.

waterfall(F,T,abs(S(:,:,1))')
helperGraphicsOpt(1)

waterfall(F,T,abs(S(:,:,2))')
helperGraphicsOpt(2)

waterfall(F,T,abs(S(:,:,3))')
helperGraphicsOpt(3)

다음 헬퍼 함수는 현재 좌표축의 모양 및 동작을 설정합니다.

function helperGraphicsOpt(ChannelId)
ax = gca;
ax.XDir = 'reverse';
ax.ZLim = [0 30];
ax.Title.String = ['Input Channel: ' num2str(ChannelId)];
ax.XLabel.String = 'Frequency (Hz)';
ax.YLabel.String = 'Time (seconds)';
ax.View = [30 45];
end

입력 인수

x — 입력 신호

벡터 | 행렬 | timetable형

입력 신호로, 벡터, 행렬 또는 MATLAB^®timetable형으로 지정됩니다.

참고

x와 s의 길이가 같도록 하려면 (length(x)-noverlap)/(length(window)-noverlap)의 값이 정수여야 합니다. Window를 사용하여 window의 길이를 지정하고 OverlapLength를 사용하여 noverlap을 지정합니다.

• 입력값에 여러 채널이 있을 경우 x를 행렬로 지정하십시오. 여기서 각 열은 채널에 해당합니다.

• 타임테이블 입력값의 경우 x는 균일한 간격으로 증가하는 유한한 행 시간값을 포함해야 합니다. 타임테이블에 누락되거나 중복된 시간 지점이 포함된 경우 누락되거나 중복되거나 불균일하게 분포된 시간값이 포함된 타임테이블 정리하기에 나와 있는 팁을 활용하여 수정할 수 있습니다.

• 다중채널 타임테이블 입력값의 경우 x를 하나의 행렬을 담은 단일 변수가 있는 타임테이블 또는 각각의 변수에 열 벡터가 담긴 복수의 변수가 있는 타임테이블로 지정합니다. 모든 변수의 정밀도가 같아야 합니다.

x의 각 채널은 윈도우 길이보다 더 크거나 같아야 합니다.

예: chirp(0:1/4e3:2,250,1,500,"quadratic")은 단일채널 처프를 지정합니다.

예: timetable(rand(5,2),SampleRate=1)은 4초 동안 1Hz로 샘플링된 2채널 확률 변수를 지정합니다.

예: timetable(rand(5,1),rand(5,1),SampleRate=1)은 4초 동안 1Hz로 샘플링된 2채널 확률 변수를 지정합니다.

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

fs — 샘플 레이트

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

샘플 레이트로, 양의 스칼라로 지정됩니다. 이 인수는 x가 벡터 또는 행렬일 때만 적용됩니다.

데이터형: double | single

ts — 샘플 시간

duration형 스칼라

샘플 시간으로, duration형 스칼라로 지정됩니다. 이 인수는 x가 벡터 또는 행렬일 때만 적용됩니다.

예: seconds(1)은 연속 신호 샘플 간의 시간 차분이 1초임을 나타내는 duration형 스칼라입니다.

데이터형: duration

이름-값 인수

선택적 인수 쌍을 Name1=Value1,...,NameN=ValueN으로 지정합니다. 여기서 Name은 인수 이름이고 Value는 대응값입니다. 이름-값 인수는 다른 인수 뒤에 와야 하지만, 인수 쌍의 순서는 상관없습니다.

예: Window=hamming(100),OverlapLength=50,FFTLength=128은 인접 세그먼트와 128개 점 FFT 사이에 50개 샘플 중첩을 가진 100-샘플 해밍 윈도우를 사용하여 데이터에 윈도우를 적용합니다.

R2021a 이전 릴리스에서는 쉼표를 사용하여 각 이름과 값을 구분하고 Name을 따옴표로 묶으십시오.

예: 'Window',hamming(100),'OverlapLength',50,'FFTLength',128은 인접 세그먼트와 128개 점 FFT 사이에 50개 샘플 중첩을 가진 100-샘플 해밍 윈도우를 사용하여 데이터에 윈도우를 적용합니다.

Window — 스펙트럼 윈도우

hann(128,"periodic") (디폴트 값) | 벡터

스펙트럼 윈도우로, 벡터로 지정됩니다. 윈도우를 지정하지 않거나 빈 값으로 지정할 경우 함수는 길이가 128인 핸 윈도우를 사용합니다. Window의 길이는 2보다 크거나 같아야 합니다.

사용 가능한 윈도우 목록은 윈도우 항목을 참조하십시오.

예: hann(N+1) 및 (1-cos(2*pi*(0:N)'/N))/2는 모두 길이가 N + 1인 핸 윈도우를 지정합니다.

데이터형: double | single

OverlapLength — 중첩 샘플의 개수

윈도우 길이의 75% (디폴트 값) | 음이 아닌 정수

중첩된 샘플의 개수로, Window 길이보다 작은 양의 정수로 지정됩니다. OverlapLength를 생략하거나 빈 값으로 지정할 경우 윈도우 길이의 75%보다 작은 가장 큰 정수로 설정됩니다. 디폴트 핸 윈도우의 경우 96개 샘플이 됩니다.

데이터형: double | single

FFTLength — DFT 지점 개수

128 (디폴트 값) | 양의 정수

DFT를 적용할 지점의 개수로, 양의 정수로 지정됩니다. 이 값이 윈도우 길이보다 크거나 같아야 합니다. 입력 신호의 길이가 DFT 길이보다 짧을 경우 데이터가 0으로 채워집니다.

데이터형: double | single

FrequencyRange — STFT 주파수 범위

"centered" (디폴트 값) | "twosided" | "onesided"

STFT 주파수 범위로, "centered", "twosided" 또는 "onesided"로 지정됩니다.

• "centered" — 중심이 맞춰진 양측 STFT를 계산합니다. FFTLength가 짝수이면 s는 구간 (–π, π] rad/sample에 대해 계산됩니다. FFTLength가 홀수이면 s는 구간 (–π, π) rad/sample에 대해 계산됩니다. 시간 정보를 지정하면 구간은 각각 (–f_s, f_s/2] cycles/unit time 및 (–f_s, f_s/2) cycles/unit time이 됩니다. 여기서 f_s는 유효한 샘플 레이트입니다.

• "twosided" — 구간 [0, 2π) rad/sample에 대해 양측 STFT를 계산합니다. 시간 정보를 지정하면 구간은 [0, f_s) cycles/unit time이 됩니다.

• "onesided" — 단측 STFT를 계산합니다. FFTLength가 짝수이면 s는 구간 [0, π] rad/sample에 대해 계산됩니다. FFTLength가 홀수이면 s는 구간 [0, π) rad/sample에 대해 계산됩니다. 시간 정보를 지정하면 구간은 각각 [0, f_s/2] cycles/unit time 및 [0, f_s/2) cycles/unit time이 됩니다. 여기서 f_s는 유효한 샘플 레이트입니다. 이 옵션은 실수 신호에만 유효합니다.

  참고

  이 인수가 "onesided"로 설정된 경우에는 stft 함수가 값을 양의 나이퀴스트 범위로 출력하고 총 전력을 보존하지 않습니다.

예제는 STFT 주파수 범위 항목을 참조하십시오.

데이터형: char | string

OutputTimeDimension — 출력 시간 차원

"acrosscolumns" (디폴트 값) | "downrows"

출력 시간 차원으로, "acrosscolumns" 또는 "downrows"로 지정됩니다. s의 시간 차원을 행 방향으로 적용하고 주파수 차원을 열 방향으로 적용하려면 이 값을 "downrows"로 설정하십시오. s의 시간 차원을 열 방향으로 적용하고 주파수 차원을 행 방향으로 적용하려면 이 값을 "acrosscolumns"로 설정하십시오. 출력 인수 없이 함수를 호출하면 이 입력값은 무시됩니다.

출력 인수

s — 단시간 푸리에 변환

행렬 | 3차원 배열

단시간 푸리에 변환으로, 행렬 또는 3차원 배열로 반환됩니다. 시간은 s의 열 방향으로 증가하고, 주파수는 행 방향으로 증가합니다. 세 번째 차원이 있을 경우 이는 입력 채널에 대응합니다.

• 신호 x에 N_x개의 시간 샘플이 있을 경우 s에는 k개의 열이 있고 k = ⌊(N_x–L)/(M–L)⌋이 성립합니다. 여기서 M은 Window의 길이이고, L은 OverlapLength이며, ⌊ ⌋ 기호는 바닥 함수를 나타냅니다.

• s의 행 개수는 FFTLength에 지정된 값과 같습니다.

데이터형: double | single

f — 주파수

벡터

STFT가 계산되는 주파수로, 벡터로 반환됩니다.

데이터형: double | single

t — 시점

벡터

시점으로, 벡터로 반환됩니다. t에는 단시간 파워 스펙트럼 추정값을 계산하는 데 사용되는 데이터 세그먼트 중심에 대응하는 시간 값이 포함됩니다.

• 샘플 레이트 fs를 사용하는 경우 벡터는 초 단위의 시간 값을 갖습니다.

• 샘플 레이트 ts를 사용하는 경우 벡터는 입력값과 동일한 시간 형식을 갖는 duration형 배열입니다.

• 시간 정보를 사용하지 않는 경우 벡터는 샘플 번호를 갖게 됩니다.

데이터형: double | single

세부 정보

단시간 푸리에 변환

단시간 푸리에 변환(STFT)은 비정상 신호의 주파수 성분이 시간의 흐름에 따라 어떻게 변하는지 분석하는 데 쓰입니다. STFT의 크기 제곱을 신호의 스펙트로그램 시간-주파수 표현이라고 합니다. 스펙트로그램에 대한 자세한 내용과 Signal Processing Toolbox™ 함수를 사용하여 스펙트로그램을 계산하는 방법은 Signal Processing Toolbox를 사용한 스펙트로그램 계산 항목을 참조하십시오.

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

$$k=\lfloor \frac{N_x-L}{M-L}\rfloor$$

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

STFT 행렬 $$X(f)=\left[\begin{array}{ccccc}{X}_1(f)& X_2(f)& X_3(f)& \cdots & X_k(f)\end{array}\right]$$의 m번째 열에는 시간 mR을 중심으로 하는 윈도우가 적용된 데이터의 DFT가 포함되어 있습니다.

$$X_m(f)={\displaystyle \sum _{n=-\infty }^{\infty }x(n)g(n-mR)e^{-j2\pi fn}}.$$

참고 문헌

[1] Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.

[2] Sharpe, Bruce. Invertibility of Overlap-Add Processing. https://gauss256.github.io/blog/cola.html, accessed July 2019.

[3] Smith, Julius Orion. Spectral Audio Signal Processing. https://ccrma.stanford.edu/~jos/sasp/, online book, 2011 edition, accessed Nov 2018.

확장 기능

tall형 배열

메모리에 담을 수 없을 정도로 많은 행을 가진 배열을 계산할 수 있습니다.

stft 함수는 tall형 배열을 지원하지만 다음과 같은 사용법 관련 참고 및 제한 사항이 있습니다.

OutputTimeDimension은 항상 지정되고 "downrows"로 설정되어야 합니다.

자세한 내용은 tall형 배열 항목을 참조하십시오.

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

GPU 코드 생성

GPU Coder™를 사용하여 NVIDIA® GPU용 CUDA® 코드를 생성할 수 있습니다.

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

코드 생성에는 타임테이블이 지원되지 않습니다.

스레드 기반 환경

MATLAB®의 backgroundPool을 사용해 백그라운드에서 코드를 실행하거나 Parallel Computing Toolbox™의 ThreadPool을 사용해 코드 실행 속도를 높일 수 있습니다.

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

• 출력 인수가 없는 구문은 지원되지 않습니다.

자세한 내용은 스레드 기반 환경에서 MATLAB 함수 실행하기 항목을 참조하십시오.

GPU 배열

Parallel Computing Toolbox™를 사용해 GPU(그래픽스 처리 장치)에서 실행하여 코드 실행 속도를 높일 수 있습니다.

이 함수는 GPU 배열을 완전히 지원합니다. 자세한 내용은 GPU에서 MATLAB 함수 실행하기 (Parallel Computing Toolbox) 항목을 참조하십시오.

버전 내역

R2019a에 개발됨

참고 항목

함수

• dlstft | istft | iscola | pspectrum | spectrogram | stftmag2sig

객체

• stftLayer

도움말 항목

• Signal Processing Toolbox를 사용한 스펙트로그램 계산
• Time-Frequency Gallery
• 시간-주파수 분석에 대한 실용적 입문