시간-주파수 분석 및 연속 웨이블릿 변환

이 예제는 연속 웨이블릿 변환의 가변 시간-주파수 분해능을 통해 선명한 시간-주파수 표현을 얻을 수 있는 방법을 보여줍니다.

CWT(연속 웨이블릿 변환)는 비정상 신호 분석에 최적인 시간-주파수 변환입니다. 비정상인 신호는 해당 신호의 주파수 영역 표현이 시간이 지남에 따라 달라진다는 것을 의미합니다. 심전도, 오디오 신호, 지진 데이터 및 기후 데이터 등 많은 신호는 비정상적(nonstationary)입니다.

쌍곡 처프 불러오기

두 개의 쌍곡 처프를 갖는 신호를 불러옵니다. 데이터는 2048Hz로 샘플링됩니다. 첫 번째 처프는 0.1초와 0.68초 사이에 활성화되며 두 번째 처프는 0.1초와 0.75초 사이에 활성화됩니다. 시간 $$t$$에 첫 번째 처프의 순시 주파수(Hz 단위)는 $$\frac{15\pi }{{(0.8-t)}^2}/2\pi$$입니다. 시간 $$t$$에 두 번째 처프의 순시 주파수는 $$\frac{5\pi }{{(0.8-t)}^2}/2\pi$$입니다. 신호를 플로팅합니다.

load hychirp
plot(t,hychirp)
grid on
title("Signal")
axis tight
xlabel("Time (s)")
ylabel('Amplitude')

시간-주파수 분석: 푸리에 변환

푸리에 변환(FT)은 신호에 존재하는 주파수 성분을 식별할 때 매우 우수합니다. 그러나 FT는 주파수 성분이 발생하는 시기를 식별하지 못합니다.

신호의 크기 스펙트럼을 플로팅합니다. 0Hz와 200Hz 사이의 영역을 확대합니다.

sigLen = numel(hychirp);
fchirp = fft(hychirp);
fr = Fs*(0:1/Fs:1-1/Fs);
plot(fr(1:sigLen/2),abs(fchirp(1:sigLen/2)),"x-")
xlabel("Frequency (Hz)")
ylabel("Amplitude")
axis tight
grid on
xlim([0 200])

시간-주파수 분석: 단시간 푸리에 변환

푸리에 변환은 시간 정보를 제공하지 않습니다. 주파수 변화가 일어나는 시기를 결정하기 위해 STFT(단시간 푸리에 변환) 접근법은 신호를 서로 다른 청크로 나누고 각 청크에 FT를 수행합니다. 시간-주파수 평면의 STFT 타일링은 여기에 표시됩니다.

STFT는 신호 이벤트가 발생하는 시간 및 주파수 두 개 모두에 대한 일부 정보를 제공합니다. 그러나 윈도우(세그먼트) 크기 선택이 중요합니다. STFT를 사용하는 시간-주파수 분석의 경우, 더 짧은 윈도우 크기를 선택하면 주파수 분해능이 낮아지는 대신 우수한 시간 분해능을 얻는 데 도움이 됩니다. 반대로 더 큰 윈도우를 선택하면 시간 분해능이 낮아지는 대신 우수한 주파수 분해능을 얻는 데 도움이 됩니다.

윈도우 크기를 선택하고 나면 해당 크기는 전체 분석이 이루어지는 동안 내내 고정된 채로 유지됩니다. 신호에서 기대하는 주파수 성분을 추정할 수 있으면 해당 정보를 사용하여 분석을 위한 윈도우 크기를 선택할 수 있습니다.

이들의 초기 시점에 두 처프의 순시 주파수는 대략 5Hz 및 15Hz입니다. 헬퍼 함수 helperPlotSpectrogram을 사용하여 시간 윈도우 크기가 200밀리초인 신호의 스펙트로그램을 플로팅합니다. helperPlotSpectrogram에 대한 소스 코드는 부록에 수록되어 있습니다. 헬퍼 함수는 스펙트로그램에 대한 순시 주파수를 검은색 파선 세그먼트로 플로팅합니다. 순시 주파수는 신호에서 조기에 확인되지만 나중에도 확인되지는 않습니다.

helperPlotSpectrogram(hychirp,t,Fs,200)

이제 helperPlotSpectrogram을 사용하여 시간 윈도우 크기가 50밀리초인 신호의 스펙트로그램을 플로팅합니다. 이제 신호의 후반에 발생하는 더 높은 주파수가 확인되지만 신호 시작 부분의 더 낮은 주파수는 확인되지 않습니다.

helperPlotSpectrogram(hychirp,t,Fs,50)

쌍곡 처프와 같은 비정상 신호의 경우, STFT 사용에 문제가 있습니다. 해당 문제는 단일 윈도우 크기라서 해당 신호의 전체 주파수 내용을 확인할 수 없다는 것입니다.

시간-주파수 분석: 연속 웨이블릿 변환

CWT(연속 웨이블릿 변환)는 STFT 고유의 분해능 문제를 극복하기 위해 만들어졌습니다. 시간-주파수 평면의 CWT 타일화가 아래에 표시되어 있습니다.

많은 실제 신호가 긴 스케일에서 발생하여 느리게 진동하는 성분을 가지고, 그러면서도 고주파수 이벤트가 갑자기 발생하거나 과도적인 경향을 가질 수 있기 때문에 이 평면의 CWT 타일화가 유용합니다. 그러나 고주파수가 오랜 시간 지속되는 것이 당연한 경우 CWT 사용은 적절하지 않을 것입니다. 시간 분해능을 얻지 않으면 주파수 분해능이 나빠질 것입니다. 그런데 그렇지 않은 경우도 꽤 자주 발생합니다. 사람의 음성 체계는 이 방식으로 작동합니다. 즉, 더 낮은 주파수에서는 주파수 국소화가 훨씬 더 좋고 높은 주파수에서는 시간 국소화가 더 좋습니다.

CWT의 스케일로그램을 플로팅합니다. 스케일로그램은 시간과 주파수의 함수로 플로팅된 CWT의 절댓값입니다. CWT의 주파수가 로그이므로 플롯은 로그 주파수 축을 사용합니다. 신호에 두 쌍곡 처프가 존재하는 것을 스케일로그램에서 명확히 확인할 수 있습니다. CWT를 사용하면 세그먼트 길이의 선택을 걱정하지 않고 신호의 지속 시간 동안 내내 순시 주파수를 정확히 추정할 수 있습니다.

cwt(hychirp,Fs)

흰색 파선은 영향 원뿔을 표시합니다. 영향 원뿔은 스케일로그램에서 경계 효과의 영향을 받을 가능성이 있는 영역을 보여줍니다. 자세한 내용은 Boundary Effects and the Cone of Influence 항목을 참조하십시오.

웨이블릿 계수의 크기가 얼마나 급격하게 증가하는지에 대한 감각을 얻으려면 헬퍼 함수 helperPlotScalogram3d를 사용하여 스케일로그램을 3차원 표면으로 플로팅합니다. helperPlotScalogram3d에 대한 소스 코드는 부록에 수록되어 있습니다.

helperPlotScalogram3d(hychirp,Fs)

헬퍼 함수 helperPlotScalogram을 사용하여 신호 및 순시 주파수의 스케일로그램을 플로팅합니다. helperPlotScalogram에 대한 소스 코드는 부록에 수록되어 있습니다. 순시 주파수는 스케일로그램 특징에 적합합니다.

helperPlotScalogram(hychirp,Fs)

부록 – 헬퍼 함수

helperPlotSpectrogram

function helperPlotSpectrogram(sig,t,Fs,timeres)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.

[px,fx,tx] = pspectrum(sig,Fs,"spectrogram",TimeResolution=timeres/1000);
hp = pcolor(tx,fx,20*log10(abs(px)));
hp.EdgeAlpha = 0;
ylims = hp.Parent.YLim;
yticks = hp.Parent.YTick;
cl = colorbar;
cl.Label.String = "Power (dB)";
axis tight
hold on
title("Time Resolution: "+num2str(timeres)+" ms")
xlabel("Time (s)")
ylabel("Hz");
dt  = 1/Fs;
idxbegin = round(0.1/dt);
idxend1 = round(0.68/dt);
idxend2 = round(0.75/dt);
instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);
instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);
plot(t(idxbegin:idxend1),(instfreq1(idxbegin:idxend1)),'k--');
hold on;
plot(t(idxbegin:idxend2),(instfreq2(idxbegin:idxend2)),'k--');
ylim(ylims);
hp.Parent.YTick = yticks;
hp.Parent.YTickLabels = yticks;
hold off
end

helperPlotScalogram

function helperPlotScalogram(sig,Fs)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.
[cfs,f] = cwt(sig,Fs);

sigLen = numel(sig);
t = (0:sigLen-1)/Fs;

hp = pcolor(t,log2(f),abs(cfs));
hp.EdgeAlpha = 0;
ylims = hp.Parent.YLim;
yticks = hp.Parent.YTick;
cl = colorbar;
cl.Label.String = "Magnitude";
axis tight
hold on
title("Scalogram and Instantaneous Frequencies")
xlabel("Seconds");
ylabel("Hz");
dt  = 1/2048;
idxbegin = round(0.1/dt);
idxend1 = round(0.68/dt);
idxend2 = round(0.75/dt);
instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);
instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);
plot(t(idxbegin:idxend1),log2(instfreq1(idxbegin:idxend1)),"k--");
hold on;
plot(t(idxbegin:idxend2),log2(instfreq2(idxbegin:idxend2)),"k--");
ylim(ylims);
hp.Parent.YTick = yticks;
hp.Parent.YTickLabels = 2.^yticks;
end

helperPlotScalogram3d

function helperPlotScalogram3d(sig,Fs)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.
figure
[cfs,f] = cwt(sig,Fs);

sigLen = numel(sig);
t = (0:sigLen-1)/Fs;
surface(t,f,abs(cfs));
xlabel("Time (s)")
ylabel("Frequency (Hz)")
zlabel("Magnitude")
title("Scalogram In 3-D")
set(gca,Yscale="log")
shading interp
view([-40 30])
end

참고 항목

앱

• 웨이블릿 시간-주파수 분석기

함수

• cwt | cwtfilterbank | waveletScattering | waveletScattering2

관련 예제

• Continuous Wavelet Analysis of Modulated Signals
• Using Wavelet Time-Frequency Analyzer App
• 웨이블릿 분석 및 딥러닝을 사용하여 시계열 분류하기
• Wavelet Time Scattering Classification of Phonocardiogram Data

세부 정보

• Practical Introduction to Time-Frequency Analysis Using the Continuous Wavelet Transform
• Boundary Effects and the Cone of Influence