시간-주파수 분석에 대한 실용적 입문

시간-주파수 분석을 사용하여 DTMF 신호에서 번호 식별하기

거의 모든 시변 신호는 정상(Stationary) 신호로 간주하기에 충분할 정도로 짧은 시간 구간으로 나눠질 수 있습니다. 시간-주파수 분석은 신호를 이러한 짧은 구간으로 분할하고 슬라이딩 윈도우에 대한 스펙트럼을 추정하는 방식으로 수행되는 것이 가장 일반적입니다. pspectrum 함수를 'spectrogram' 옵션과 함께 사용하면 각 슬라이딩 윈도우에 대해 FFT 기반 스펙트럼 추정값을 계산하고 신호의 주파수 성분이 시간이 지남에 따라 어떻게 변화하는지 시각화할 수 있습니다.

디지털 전화 다이얼의 신호 시스템을 살펴보겠습니다. 이러한 시스템에서 생성되는 신호를 이중톤 다중 주파수(DTMF) 신호라고 합니다. 각각의 눌린 번호에서 생성되는 소리는 2개의 상호 배타적인 주파수 집단에서 가져온 주파수를 갖는 2개의 정현파 $-$, 즉 2개의 톤 $-$의 합계로 구성됩니다. 각각의 톤 쌍은 저주파수 무리의 주파수 하나(697Hz, 770Hz, 852Hz 또는 941Hz)와 고주파수 무리의 주파수 하나(1209Hz, 1336Hz 또는 1477Hz)를 포함하며, 고유한 기호를 나타냅니다. 다음은 전화 패드의 버튼에 할당된 주파수입니다.

DTMF 신호를 생성하고 이를 들어봅니다.

[tones, Fs] = helperDTMFToneGenerator();
p = audioplayer(tones,Fs,16);
play(p)

신호를 들어보면 3자리 번호를 눌렀음을 알 수 있습니다. 그러나 어떤 번호인지는 알 수 없습니다. 다음으로, 신호를 시간 영역과 650Hz~1500Hz 대역의 주파수 영역에서 시각화합니다. 사각 윈도우를 사용하고 주파수 분해능을 개선하기 위해 pspectrum 함수의 'Leakage' 파라미터를 1로 설정합니다.

N = numel(tones);
t = (0:N-1)/Fs; 
subplot(2,1,1)
plot(1e3*t,tones)
xlabel('Time (ms)')
ylabel('Amplitude')
title('DTMF Signal')
subplot(2,1,2)
pspectrum(tones,Fs,'Leakage',1,'FrequencyLimits',[650, 1500])

신호의 시간 영역 플롯을 보면 3개의 눌러진 버튼에 해당하는 3개의 에너지 버스트가 있음을 알 수 있습니다. 버스트 길이를 측정하려면 RMS 포락선의 펄스 폭을 구하면 됩니다.

env = envelope(tones,80,'rms');
pulsewidth(env,Fs)

ans = 3×1

    0.1041
    0.1042
    0.1047

title('Pulse Width of RMS Envelope')

길이가 각각 대략 100밀리초인 펄스 3개를 볼 수 있습니다. 그러나, 어떤 번호를 눌렀는지는 알 수 없습니다. 주파수 영역 플롯은 신호에 존재하는 주파수를 보여주기 때문에 이를 파악하는 데 도움이 됩니다.

4개의 서로 다른 주파수 대역에서 평균 주파수를 추정하여 주파수 피크를 찾습니다.

f = [meanfreq(tones,Fs,[700 800]), ...
     meanfreq(tones,Fs,[800 900]), ...
     meanfreq(tones,Fs,[900 1000]), ...
     meanfreq(tones,Fs,[1300 1400])];
round(f)

ans = 1×4

         770         852         941        1336

추정된 주파수를 전화 패드 도식에 맞춰보면, '5', '8', '0'이 눌러진 버튼임을 알 수 있습니다. 하지만, 주파수 영역 플롯은 버튼이 눌러진 순서를 파악할 수 있는 그 어떤 유형의 시간 정보도 제공하지 않습니다. 이에 대한 가능한 조합은 '580', '508', '805', '850', '085' 또는 '058'입니다. 이 퍼즐을 풀기 위해 pspectrum 함수를 사용하여 스펙트로그램을 계산하고 신호의 주파수 성분이 시간에 따라 어떻게 달라지는지 관찰합니다.

기본 주파수 성분만 시각화하기 위해, 650Hz~1500Hz 대역에서 스펙트로그램을 계산하고 전력 수준이 $-$10dB 아래인 성분을 제거합니다. 톤의 지속 시간과 시간 영역에서의 해당 위치를 보기 위해 0% 중첩을 사용합니다.

pspectrum(tones,Fs,'spectrogram','Leakage',1,'OverlapPercent',0, ...
    'MinThreshold',-10,'FrequencyLimits',[650, 1500]);

스펙트로그램의 색은 주파수 전력 수준을 인코딩합니다. 노란색은 전력이 높은 주파수 성분을 나타내고, 파란색은 전력이 매우 낮은 주파수 성분을 나타냅니다. 진한 노란색 가로줄은 특정 주파수에서 톤이 존재함을 나타냅니다. 이 플롯은 눌러진 3개의 번호에 모두 1336Hz 톤이 존재한다는 것을 명확히 보여주므로 번호가 모두 키패드의 두 번째 열에 있다는 것을 알 수 있습니다. 플롯을 보면 가장 낮은 주파수 770Hz를 가장 먼저 눌렀음을 알 수 있습니다. 가장 높은 주파수 941Hz가 그 다음이고, 중간 주파수 852Hz가 마지막입니다. 따라서, 누른 번호는 508입니다.

시간 분해능과 주파수 분해능을 절충하여 신호에 가장 적합한 표현 구하기

pspectrum 함수는 신호를 세그먼트로 나눕니다. 세그먼트가 길수록 주파수 분해능이 향상되고, 세그먼트가 짧을수록 시간 분해능이 향상됩니다. 세그먼트 길이는 'FrequencyResolution' 파라미터와 'TimeResolution' 파라미터를 사용하여 제어할 수 있습니다. 지정된 주파수 분해능 값이나 시간 분해능 값이 없는 경우, pspectrum은 입력 신호 길이를 기준으로 시간 분해능과 주파수 분해능 간에 적절한 균형을 찾으려고 합니다.

다음에서 살펴볼 신호는 태평양 흰긴수염고래가 부르는 노래 중 짧게 반복되는 소리로 구성되었으며, 4kHz로 샘플링되었습니다.

load whaleTrill
p = audioplayer(whaleTrill,Fs,16);
play(p)

짧게 반복되는 소리 신호는 일련의 음조 펄스로 이루어집니다. 이제 시간 신호와 함께, 분해능을 지정하지 않았을 때와 시간 분해능을 10밀리초로 설정했을 때 pspectrum으로 구한 스펙트로그램을 살펴보겠습니다. 'Leakage' 파라미터를 1로 설정하여 사각 윈도우를 사용합니다. 펄스의 시간 위치를 국한하고자 하기 때문에 중첩 비율을 0으로 설정합니다. 마지막으로, $-$60dB의 'MinThreshold'를 사용하여 스펙트로그램 보기에서 배경 잡음을 제거합니다.

t = (0:length(whaleTrill)-1)/Fs;
figure
ax1 = subplot(3,1,1);
plot(t,whaleTrill)
ax2 = subplot(3,1,2);
pspectrum(whaleTrill,Fs,'spectrogram','OverlapPercent',0, ...
    'Leakage',1,'MinThreshold',-60)
colorbar(ax2,'off')
ax3 = subplot(3,1,3);
pspectrum(whaleTrill,Fs,'spectrogram','OverlapPercent',0, ...
    'Leakage',1,'MinThreshold',-60,'TimeResolution', 10e-3)
colorbar(ax3,'off')
linkaxes([ax1,ax2,ax3],'x')

pspectrum에서 선택한 47밀리초 시간 분해능은 스펙트로그램에서 짧게 반복되는 소리 펄스를 모두 국한할 만큼 크기가 작지 않습니다. 반면, 10밀리초 시간 분해능은 짧게 반복되는 각 소리 펄스를 시간 영역에서 특정 시점으로 충분히 국한할 수 있습니다. 이는 몇몇 펄스를 확대해 보면 더 명확하게 알 수 있습니다.

xlim([0.3 0.55])

이제 큰갈색박쥐(Eptesicus fuscus)가 발사한 반향정위 펄스로 구성된 신호를 불러오겠습니다. 이 신호는 7마이크로초의 샘플링 간격으로 측정되었습니다. 신호의 스펙트로그램을 분석합니다.

load batsignal
Fs = 1/DT;
figure
pspectrum(batsignal,Fs,'spectrogram')

디폴트 파라미터 값으로 구성된 스펙트로그램에는 4개의 성긴 시간-주파수 리지가 표시됩니다. 각 리지의 주파수 변동에 대한 자세한 정보를 얻기 위해 주파수 분해능 값을 3kHz로 줄입니다.

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3)

이제 주파수 리지가 주파수 영역에서 특정 주파수로 더 잘 국한되었음을 알 수 있습니다. 그러나 주파수 분해능은 시간 분해능에 반비례하기 때문에, 스펙트로그램의 시간 분해능이 훨씬 더 작습니다. 시간 윈도우를 평활화하기 위해 중첩을 99%로 설정합니다. 'MinThreshold' $-$60dB을 사용하여 원치 않는 배경 성분을 제거합니다.

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3, ...
    'OverlapPercent',99,'MinTHreshold',-60)

이 새로운 설정은 반향정위 신호의 주파수 리지 4개를 명확하게 표시하는 스펙트로그램을 생성합니다.

시간-주파수 재할당

4개의 주파수 리지를 식별할 수는 있지만, 여전히 각 리지가 여러 인접 주파수 Bin에 걸쳐 분산되어 있음을 알 수 있습니다. 이는 시간과 주파수 모두에 사용되는 윈도우 적용 방법의 누설로 인한 것입니다.

pspectrum 함수는 시간과 주파수 모두에서 각 스펙트럼 추정값에 대한 에너지 중심을 추정할 수 있습니다. 각 추정값의 에너지를 새로운 시간 및 주파수 중심에 가장 가까운 Bin에 재할당하면 윈도우의 누설을 일부 해결할 수 있습니다. 'Reassign' 파라미터를 사용하면 가능합니다. 이 파라미터를 true로 설정하면 신호의 재할당된 스펙트로그램이 계산됩니다.

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3, ...
    'OverlapPercent',99,'MinTHreshold',-60,'Reassign',true)

이제 주파수 리지가 훨씬 더 매끄러워졌으며, 시간 영역에서 특정 시간으로 더 잘 국한되었습니다. 또한 함수 fsst를 사용하여 신호 에너지를 국한할 수도 있습니다. 이에 대해서는 다음 섹션에서 설명합니다.

시간-주파수 리지(Ridge) 재생성하기

시간이 지남에 따라 주파수가 감소하는 처프 신호(Chirp Signal)와 마지막에 무엇인가 세게 부딪히는 "철퍼덕" 소리로 구성된 다음 녹음을 살펴보겠습니다.

load splat
p = audioplayer(y,Fs,16);
play(p)
pspectrum(y,Fs,'spectrogram')

시간-주파수 평면에서 리지를 추출하여 "철퍼덕" 소리의 일부를 재구성해 보겠습니다. 여기서는 fsst를 사용하여 잡음이 많은 "철퍼덕" 신호에 대한 스펙트럼을 매끄럽게 하고, tfridge를 사용하여 처프 소리의 리지를 식별하며, ifsst를 사용하여 처프를 재생성합니다. 그 과정에서, 재구성된 신호의 잡음이 제거됩니다.

가우스 잡음을 "철퍼덕" 소리의 처프 부분에 추가합니다. 추가된 잡음 덕분에 저가의 마이크로 녹음한 듯한 결과를 얻게 됩니다. 시간-주파수 스펙트럼 성분을 검토합니다.

rng('default')
t = (0:length(y)-1)/Fs;
yNoise = y + 0.1*randn(size(y));
yChirp = yNoise(t<0.35);
pspectrum(yChirp,Fs,'spectrogram','MinThreshold',-70)

푸리에 싱크로스퀴즈 변환(Fourier Synchrosqueezed Transform), 즉 fsst를 사용하여 스펙트럼을 매끄럽게 만듭니다. fsst는 고정 시간 동안 주파수 상에서 에너지를 재할당하여 시간-주파수 평면의 특정 영역에 에너지를 국한합니다. 잡음이 있는 처프의 싱크로스퀴즈 변환을 계산하고 플로팅합니다.

fsst(yChirp,Fs,'yaxis')

처프가 시간-주파수 평면에서 국한된 리지로 나타납니다. tfridge를 사용하여 리지를 식별합니다. 변환과 함께 리지를 플로팅합니다.

[sst,f] = fsst(yChirp,Fs); 
[fridge, iridge] = tfridge(sst,f,10);
helperPlotRidge(yChirp,Fs,fridge);

다음으로, 리지 인덱스 벡터 iridge를 사용하여 처프 신호를 재생성합니다. 리지의 양쪽에 하나의 Bin을 포함시킵니다. 재생성된 신호의 스펙트로그램을 플로팅합니다.

yrec = ifsst(sst,kaiser(256,10),iridge,'NumFrequencyBins',1);
pspectrum(yrec,Fs,'spectrogram','MinThreshold',-70)

리지를 재생성함으로써 신호에서 잡음이 제거되었습니다. 잡음이 있는 신호와 잡음이 제거된 신호를 연속으로 재생하여 차이를 들어보십시오.

p = audioplayer([yChirp;zeros(size(yChirp));yrec],Fs,16);
play(p);

전력 측정하기

일반적인 레이다 파형인 복소 선형 주파수 변조(LFM) 펄스를 살펴보겠습니다. 1.27마이크로초의 시간 분해능과 90% 중첩을 사용하여 신호의 스펙트로그램을 계산합니다.

Fs = 1e8;
bw = 60e6;
t = 0:1/Fs:10e-6;
IComp = chirp(t,-bw/2,t(end), bw/2,'linear',90)+0.15*randn(size(t));
QComp = chirp(t,-bw/2,t(end), bw/2,'linear',0) +0.15*randn(size(t));
IQData = IComp + 1i*QComp;

segmentLength = 128;
pspectrum(IQData,Fs,'spectrogram','TimeResolution',1.27e-6,'OverlapPercent',90)

스펙트로그램 계산에 사용된 파라미터는 LFM 신호에 대한 명확한 시간-주파수 표현을 제공합니다. pspectrum은 파워 스펙트로그램을 계산하는데, 이는 색 값이 실제 전력 수준(단위: dB)에 대응함을 의미합니다. 컬러바를 보면 신호의 전력 수준이 약 $-$4dB임을 알 수 있습니다.

로그 주파수 스케일 시각화

특정 응용 사례에서는 로그 주파수 스케일을 사용해 신호의 스펙트로그램을 시각화하는 것이 더 바람직할 수 있습니다. y축의 YScale 속성을 변경하면 그렇게 할 수 있습니다. 예를 들어, 1kHz로 샘플링된 로그 처프(Chirp)가 있다고 가정하겠습니다. 처프의 주파수는 10초 동안 10Hz에서 400Hz로 증가합니다.

Fs = 1e3;                    
t = 0:1/Fs:10;               
fo = 10;                     
f1 = 400;                    
y = chirp(t,fo,10,f1,'logarithmic');
pspectrum(y,Fs,'spectrogram','FrequencyResolution',1, ...
    'OverlapPercent',90,'Leakage',0.85,'FrequencyLimits',[1 Fs/2])

처프의 스펙트로그램은 주파수 스케일이 로그일 때 직선으로 나타납니다.

ax = gca;
ax.YScale = 'log';

3차원 폭포 플롯 시각화

view 명령을 사용하여 신호의 스펙트로그램을 3차원 폭포 플롯으로 시각화할 수 있습니다. 또한 colormap 함수를 사용하여 표시 색을 변경할 수도 있습니다.

Fs = 10e3;
t = 0:1/Fs:2;
x1 = vco(sawtooth(2*pi*t,0.5),[0.1 0.4]*Fs,Fs);
pspectrum(x1,Fs,'spectrogram','Leakage',0.8)

view(-45,65)
colormap bone

지속성 스펙트럼을 사용하여 간섭 찾기

신호의 지속성 스펙트럼은 주어진 주파수가 신호에 존재하는 시간 비율을 보여주는 시간-주파수 보기입니다. 지속성 스펙트럼은 전력-주파수 영역의 히스토그램입니다. 신호가 전개될 때 신호의 특정 주파수가 오래 지속될수록 그 주파수가 차지하는 시간 비율이 높아지기 때문에 해당 주파수 색이 디스플레이에 더 밝거나 "뜨겁게" 표시됩니다. 지속성 스펙트럼을 사용하여 다른 신호에 숨어 있는 신호를 식별할 수 있습니다.

광대역 신호에 포함된 협대역 간섭 신호를 살펴보겠습니다. 500초 동안 1kHz로 샘플링된 처프를 생성합니다. 측정하는 동안 처프의 주파수는 180Hz에서 220Hz로 증가합니다.

fs = 1000;
t = (0:1/fs:500)';

x = chirp(t,180,t(end),220) + 0.15*randn(size(t));

또한 신호에는 진폭이 0.05인 210Hz 간섭도 포함되는데, 이 간섭은 총 신호 지속 시간의 1/6 동안에만 존재합니다.

idx = floor(length(x)/6);
x(1:idx) = x(1:idx) + 0.05*cos(2*pi*t(1:idx)*210);

100Hz~290Hz 구간에 대한 신호의 파워 스펙트럼을 계산합니다. 약한 정현파가 처프에 의해 가려집니다.

pspectrum(x,fs,'FrequencyLimits',[100 290])

신호의 지속성 스펙트럼을 계산합니다. 이제 두 신호 성분이 명확히 보입니다.

figure
colormap parula
pspectrum(x,fs,'persistence','FrequencyLimits',[100 290],'TimeResolution',1)

결론

이 예제에서는 pspectrum 함수를 사용하여 시간-주파수 분석을 수행하는 방법과 스펙트로그램 데이터와 전력 수준을 해석하는 방법에 대해 배웠습니다. 그리고 신호에 대한 이해를 돕기 위해 시간 분해능과 주파수 분해능을 변경하는 방법과 fsst, ifsst, tfridge를 사용하여 스펙트럼을 매끄럽게 만들고 시간-주파수 리지를 추출하는 방법을 배웠습니다. 또한 스펙트로그램 플롯을 구성하여 로그 주파수 스케일과 3차원 시각화를 구현하는 방법도 배웠습니다. 마지막으로, 지속성 스펙트럼을 계산하여 간섭 신호를 찾는 방법을 알아보았습니다.

부록

이 예제에서는 다음 헬퍼 함수가 사용되었습니다.