Main Content

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

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

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

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

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

2초 동안 1kHz로 샘플링된 2차 처프(Chirp)를 생성합니다. 순시 주파수는 t=0초에서 100Hz, 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);

Figure contains an axes object. The axes object with title Short-Time Fourier Transform, xlabel Time (s), ylabel Frequency (Hz) contains an object of type image.

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

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

stft 디폴트 값

spectrogramstft 함수를 사용하여 신호의 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=Nx/4.5인 세그먼트로 나눕니다. 여기서 Nx는 신호의 길이입니다. 해밍 윈도우로 각 세그먼트에 윈도우를 적용합니다.

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

  • max(256,2log2M)개의 점을 사용하여 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

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)

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

신호의 단측, 양측 및 중심이 맞춰진 단시간 푸리에 변환을 계산합니다. 모든 경우에, 형태 인자 β=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

Figure contains 3 axes objects. Axes object 1 with title 'onesided': [0.000, 2.500] kHz, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image. Axes object 2 with title 'twosided': [0.000, 4.975] kHz, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image. Axes object 3 with title 'centered': [-2.475, 2.500] kHz, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image.

계산을 반복하되, 이번에는 카이저 윈도우의 길이를 홀수인 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

Figure contains 3 axes objects. Axes object 1 with title 'onesided': [0.000, 2.488] kHz, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image. Axes object 2 with title 'twosided': [0.000, 4.975] kHz, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image. Axes object 3 with title 'centered': [-2.488, 2.488] kHz, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image.

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)

Figure contains an axes object. The axes object with title Input Channel: 1, xlabel Frequency (Hz), ylabel Time (seconds) contains an object of type patch.

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

Figure contains an axes object. The axes object with title Input Channel: 2, xlabel Frequency (Hz), ylabel Time (seconds) contains an object of type patch.

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

Figure contains an axes object. The axes object with title Input Channel: 3, xlabel Frequency (Hz), ylabel Time (seconds) contains an object of type patch.

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

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

입력 인수

모두 축소

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

참고

xs의 길이가 같도록 하려면 (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
복소수 지원 여부:

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

데이터형: double | single

샘플 시간으로, 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-샘플 해밍 윈도우를 사용하여 데이터에 윈도우를 적용합니다.

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

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

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

데이터형: double | single

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

데이터형: double | single

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

데이터형: double | single

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

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

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

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

    참고

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

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

데이터형: char | string

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

출력 인수

모두 축소

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

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

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

데이터형: double | single

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

데이터형: double | single

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

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

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

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

데이터형: double | single

세부 정보

모두 축소

단시간 푸리에 변환

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

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

k=NxLML

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

STFT 행렬 X(f)=[X1(f)X2(f)X3(f)Xk(f)]의 m번째 열에는 시간 mR을 중심으로 하는 윈도우가 적용된 데이터의 DFT가 포함되어 있습니다.

Xm(f)=n=x(n)g(nmR)ej2πfn.

  • 단시간 푸리에 변환은 가역적입니다. 역변환 프로세스는 윈도우가 적용된 세그먼트를 중첩 가산(overlap-add)하여 윈도우 모서리에서 신호 감쇠를 보정합니다. 자세한 내용은 Inverse Short-Time Fourier Transform 항목을 참조하십시오.

  • istft 함수는 신호의 STFT를 역변환합니다.

  • 특정한 상황에서는 신호의 "완벽한 복원"이 가능할 수 있습니다. 자세한 내용은 Perfect Reconstruction 항목을 참조하십시오.

  • stftmag2sig 함수는 STFT의 크기로부터 복원된 신호의 추정값을 반환합니다.

완전한 복원

일반적으로 입력 신호의 STFT를 계산하고 이에 대한 역변환을 계산하더라도 원래의 값이 완전히 복원되지는 않습니다. ISTFT 출력이 원본 입력 신호와 최대한 일치되도록 하려면 신호와 윈도우는 다음 조건을 충족해야 합니다.

  • 입력 크기 — istft를 사용하여 stft의 출력값에 대한 역변환을 계산했을 때 그 결과가 입력 신호 x와 동일한 길이가 되게 하려면

    k = NxLML

    의 값이 정수여야 합니다. 수식에서 Nx는 신호의 길이이고, M은 윈도우의 길이이며, L은 중첩 길이입니다.

  • COLA 준수 — 신호의 단시간 푸리에 변환을 수정하지 않았다면 COLA 준수 윈도우를 사용해야 합니다.

  • 채우기 — 입력 신호의 길이로 인해 k의 값이 정수가 되지 않도록 만드는 경우 단시간 푸리에 변환의 계산에 앞서 신호를 0으로 채워야 합니다. 신호에 대한 역변환을 계산한 다음 추가된 0을 제거해야 합니다.

참고 문헌

[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.

확장 기능

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

버전 내역

R2019a에 개발됨