cwt

디폴트 값을 사용하여 음성 샘플의 연속 웨이블릿 변환을 구합니다.

load mtlb;
w = cwt(mtlb);

음성 샘플을 불러옵니다.

load mtlb

파일 mtlb.mat를 불러오면 음성 신호 mtlb와 샘플 레이트 Fs가 작업 공간으로 들어옵니다. 범프 웨이블릿을 사용하여 얻은 음성 샘플의 스케일로그램을 표시합니다.

load mtlb
cwt(mtlb,"bump",Fs)

디폴트 Morse 웨이블릿을 사용하여 얻은 스케일로그램과 비교합니다.

cwt(mtlb,Fs)

고베 지진 데이터의 CWT를 구합니다. 이 데이터는 1995년 1월 16일 20:56:51(GMT)부터 51분 동안 호주 Tasmania University, Hobart에서 기록한 지진계(수직 가속도, nm/sq.sec) 측정값입니다. 샘플링 주파수는 1Hz입니다.

load kobe

지진 데이터를 플로팅합니다.

plot((1:numel(kobe))./60,kobe)
xlabel("Time (mins)")
ylabel("Vertical Acceleration (nm/s^2)")
title("Kobe Earthquake Data")
grid on
axis tight

CWT, 주파수, 영향 원뿔을 구합니다.

[wt,f,coi] = cwt(kobe,1);

영향 원뿔을 포함하여 스케일로그램을 확인합니다.

cwt(kobe,1)

샘플링 주파수 대신 샘플링 주기를 지정하여 CWT, 시간 주기 및 영향 원뿔을 구합니다.

[wt,periods,coi] = cwt(kobe,minutes(1/60));

샘플링 주기를 지정할 때 생성된 스케일로그램을 확인합니다.

cwt(kobe,minutes(1/60))

32Hz 및 64Hz의 주파수를 갖고 진폭이 서로 다른 두 개의 복소수 지수를 생성합니다. 데이터는 1000Hz로 샘플링됩니다. 두 복소수 지수는 시간적으로 서로 겹치지 않는 지지 범위를 갖습니다.

Fs = 1e3;
t = 0:1/Fs:1;
z = exp(1i*2*pi*32*t).*(t>=0.1 & t<0.3)+2*exp(-1i*2*pi*64*t).*(t>0.7);

표준편차가 0.05인 복소수 백색 가우스 잡음을 추가합니다.

wgnNoise = 0.05/sqrt(2)*randn(size(t))+1i*0.05/sqrt(2)*randn(size(t));
z = z+wgnNoise;

Morse 웨이블릿을 사용하여 cwt를 구하고 플로팅합니다.

cwt(z,Fs)

컬러바에 있는 복소수 지수 성분의 크기는 스케일이 다르더라도 본질적으로 진폭입니다. 이는 L1 정규화의 직접적인 결과입니다. 스크립트를 실행하고 데이터 커서로 각 서브플롯을 탐색하여 이를 확인할 수 있습니다.

이 예제에서는 신호에서 진동 성분의 진폭이 그에 대응하는 웨이블릿 계수의 진폭과 일치함을 보여줍니다.

시간적으로 서로 겹치지 않는 지지 범위를 갖는 두 개의 정현파로 구성된 신호를 생성합니다. 한 정현파는 주파수가 32Hz이고 진폭이 1입니다. 다른 정현파는 주파수가 64Hz이고 진폭이 2입니다. 이 신호는 1000Hz에서 1초 동안 샘플링됩니다. 신호를 플로팅합니다.

frq1 = 32;
amp1 = 1;
frq2 = 64;
amp2 = 2;

Fs = 1e3;
t = 0:1/Fs:1;
x = amp1*sin(2*pi*frq1*t).*(t>=0.1 & t<0.3)+... 
    amp2*sin(2*pi*frq2*t).*(t>0.6 & t<0.9);

plot(t,x)
grid on
xlabel("Time (sec)")
ylabel("Amplitude")
title("Signal")

신호에 적용할 수 있는 CWT 필터 뱅크를 만듭니다. 신호 성분 주파수를 알고 있으므로 필터 뱅크의 주파수 제한을 알려진 주파수가 포함된 좁은 범위로 설정합니다. 범위를 확인하기 위해 필터 뱅크에 대한 크기 주파수 응답을 플로팅합니다.

fb = cwtfilterbank(SignalLength=numel(x),SamplingFrequency=Fs,...
    FrequencyLimits=[20 100]);
freqz(fb)

cwt 및 필터 뱅크를 사용하여 신호의 스케일로그램을 플로팅합니다.

cwt(x,FilterBank=fb)

데이터 커서를 사용하여 웨이블릿 계수의 진폭이 본질적으로 정현파 성분의 진폭과 같은지 확인합니다. 결과는 다음 그림과 유사해야 합니다.

이 예제에서는 여러 시계열에 CWT를 수행할 때 CWT 필터 뱅크를 사용하여 계산 효율성을 개선할 수 있는 방법을 보여줍니다.

100×1024 행렬 x를 생성합니다. 1024개의 샘플이 있는 신호에 적합한 CWT 필터 뱅크를 생성합니다.

x = randn(100,1024);
fb = cwtfilterbank;

디폴트 설정의 cwt를 사용하여 1024개의 샘플이 있는 신호의 CWT를 구합니다. 각각 1024개의 샘플이 있는 신호 100개의 CWT 계수를 포함할 수 있는 3차원 배열을 생성합니다.

cfs = cwt(x(1,:));
res = zeros(100,size(cfs,1),size(cfs,2));

cwt 함수를 사용하여 행렬 x의 각 행에 CWT를 수행합니다. 경과 시간을 표시합니다.

tic
for k=1:100
    res(k,:,:) = cwt(x(k,:));
end
toc

Elapsed time is 0.928160 seconds.

이제 필터 뱅크의 wt 객체 함수를 사용하여 x의 각 행에 CWT를 수행합니다. 경과 시간을 표시합니다.

tic
for k=1:100
    res(k,:,:) = wt(fb,x(k,:));
end
toc

Elapsed time is 0.393524 seconds.

이 예제에서는 생성된 CUDA® 코드를 사용하여 CWT(연속 웨이블릿 변환)를 수행하기 위해 MEX 파일을 생성하는 방법을 보여줍니다.

먼저 CUDA 지원 GPU와 NVCC 컴파일러가 있는지 확인하십시오. 구성이 적절한지 확인하려면 The GPU Environment Check and Setup App (GPU Coder) 항목을 참조하십시오.

GPU 코더 구성 객체를 만듭니다.

cfg = coder.gpuConfig("mex");

1,000Hz에서 100,000개 샘플의 신호를 생성합니다. 신호는 시간적으로 서로 겹치지 않는 지지 범위를 갖는 두 개의 코사인파로 구성됩니다.

t = 0:.001:(1e5*0.001)-0.001;
x = cos(2*pi*32*t).*(t > 10 & t<=50)+ ...
    cos(2*pi*64*t).*(t >= 60 & t < 90)+ ...
    0.2*randn(size(t));

단정밀도를 사용하도록 신호를 형변환합니다. GPU 계산은 단정밀도에서 더 효율적인 경우가 많습니다. 그러나 NVIDIA® GPU가 배정밀도를 지원하는 경우 배정밀도 코드를 생성할 수도 있습니다.

x = single(x);

GPU MEX 파일과 코드 생성 리포트를 생성합니다. MEX 파일 생성을 허용하려면 다음 입력 파라미터 3개의 속성(클래스, 크기, 실수/복소수 여부)을 지정해야 합니다.

sig = coder.typeof(single(0),[1 1e5]);
wav = coder.typeof('c',[1 inf]);
sfrq = coder.typeof(0);
codegen cwt -config cfg -args {sig,wav,sfrq} -report

Code generation successful: View report

-report 플래그는 선택 사항입니다. -report를 사용하면 코드 생성 리포트가 생성됩니다. 리포트의 요약 탭에서 GPU 코드 메트릭 링크를 찾을 수 있습니다. 이 링크에서는 생성된 CUDA 커널의 개수 및 할당된 메모리 양과 같은 자세한 정보를 제공합니다.

데이터에 대해 MEX 파일을 실행하고 스케일로그램을 플로팅합니다. 플롯이 서로 겹치지 않는 두 코사인파와 일치하는지 확인합니다.

[cfs,f] = cwt_mex(x,'morse',1e3);
image("XData",t,"YData",f,"CData",abs(cfs),"CDataMapping","scaled")
set(gca,"YScale","log")
axis tight
xlabel("Time (Seconds)")
ylabel("Frequency (Hz)")
title("Scalogram of Two-Tone Signal")

_mex를 추가하지 않고 위의 CWT 명령을 실행합니다. MATLAB® 스케일로그램과 GPU MEX 스케일로그램이 동일한지 확인합니다.

[cfs2,f2] = cwt(x,'morse',1e3);
max(abs(cfs2(:)-cfs(:)))

ans = single
    7.3380e-07

이 예제에서는 출력 인수 없이 플롯을 얻을 때 CWT의 디폴트 주파수 축 레이블을 변경하는 방법을 보여줍니다.

32Hz 및 64Hz 주파수를 갖는 두 개의 사인파를 생성합니다. 데이터는 1000Hz로 샘플링됩니다. 두 사인파는 시간적으로 서로 겹치지 않는 지지 범위를 갖습니다. 표준편차가 0.05인 백색 가우스 잡음을 추가합니다. 디폴트 Morse 웨이블릿을 사용하여 CWT를 얻고 플로팅합니다.

Fs = 1e3;
t = 0:1/Fs:1;
x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
cwt(x,1000)

CWT의 주파수가 로그이므로 플롯은 로그 주파수 축을 사용합니다. MATLAB에서 로그 좌표축의 단위는 10의 거듭제곱(디케이드)입니다. cwtfreqbounds를 사용하여, 주어진 신호 길이, 샘플링 주파수, 웨이블릿에 대한 최소 및 최대 웨이블릿 대역통과 주파수를 확인할 수 있습니다.

[minf,maxf] = cwtfreqbounds(numel(x),1000);

기본적으로 주파수 눈금은 최소 주파수와 최대 주파수 사이의 10의 거듭제곱이기 때문에 10과 100에 배치되었음을 알 수 있습니다. 주파수 축 눈금을 추가하려면 다음을 사용하여 최소 주파수와 최대 주파수 사이의 로그 간격 주파수 집합을 얻으면 됩니다.

numfreq = 10;
freq = logspace(log10(minf),log10(maxf),numfreq);

그런 다음, 다음과 같은 방법으로 현재 좌표축에 대한 핸들을 가져오고 주파수 축 눈금과 레이블을 바꿉니다.

AX = gca;
AX.YTickLabelMode = "auto";
AX.YTick = freq;

CWT에서 주파수는 2의 거듭제곱으로 계산됩니다. 주파수 눈금과 눈금 레이블을 2의 거듭제곱 단위로 만들려면 다음을 수행하면 됩니다.

newplot
cwt(x,1000)
AX = gca;
freq = 2.^(round(log2(minf)):round(log2(maxf)));
AX.YTickLabelMode = "auto";
AX.YTick = freq;

이 예제에서는 플로팅을 위해 각 수준의 최대 절댓값으로 스케일로그램 값을 스케일링하는 방법을 보여줍니다.

신호를 불러오고 디폴트 스케일로그램을 표시합니다. 컬러맵을 pink(240)으로 변경합니다.

load noisdopp
cwt(noisdopp)
colormap(pink(240))

신호에 CWT를 수행하고 웨이블릿 계수와 주파수를 구합니다.

[cfs,frq] = cwt(noisdopp);

각 주파수(수준)에서 계수의 최댓값을 효율적으로 구하려면 먼저 계수의 절댓값을 전치합니다. 모든 수준에서 최솟값을 구합니다. 각 수준에서 해당 수준의 최솟값을 뺍니다.

tmp1 = abs(cfs);
t1 = size(tmp1,2);
tmp1 = tmp1';
minv = min(tmp1);
tmp1 = (tmp1-minv(ones(1,t1),:));

tmp1의 모든 수준에서 최댓값을 구합니다. 각 수준에 대해 모든 값을 해당 수준의 최댓값으로 나눕니다. 그 결과에 컬러맵의 색 개수를 곱합니다. 0인 요소를 모두 1로 설정합니다. 결과를 전치합니다.

maxv = max(tmp1);
maxvArray = maxv(ones(1,t1),:);
indx = maxvArray<eps;
tmp1 = 240*(tmp1./maxvArray);
tmp2 = 1+fix(tmp1);
tmp2(indx) = 1;
tmp2 = tmp2';

결과를 표시합니다. 스케일로그램 값은 이제 각 수준의 최대 절댓값으로 스케일링됩니다. 주파수는 선형 스케일로 표시됩니다.

t = 0:length(noisdopp)-1;
pcolor(t,frq,tmp2)
shading interp
xlabel("Time (Samples)")
ylabel("Normalized Frequency (cycles/sample)")
title("Scalogram Scaled By Level")
colormap(pink(240))
colorbar

이 예제에서는 Morse 웨이블릿의 시간과 대역폭의 곱 $$P^2$$을 증가시키면 포락선 아래에 더 많은 진동을 갖는 웨이블릿이 생성된다는 것을 알 수 있습니다. $$P^2$$이 증가하면 주파수 영역에서 웨이블릿이 좁아집니다.

두 개의 필터 뱅크를 생성합니다. 첫 번째 필터 뱅크의 디폴트 TimeBandwidth 값은 60입니다. 두 번째 필터 뱅크의 TimeBandwidth 값은 10입니다. 두 필터 뱅크의 SignalLength는 4096개 샘플입니다.

sigLen = 4096;
fb60 = cwtfilterbank(SignalLength=sigLen);
fb10 = cwtfilterbank(SignalLength=sigLen,TimeBandwidth=10);

필터 뱅크에 대한 시간 영역 웨이블릿을 구합니다.

[psi60,t] = wavelets(fb60);
[psi10,~] = wavelets(fb10);

scales 함수를 사용하여 각 필터 뱅크에 대한 모 웨이블릿(mother wavelet)을 구합니다.

sca60 = scales(fb60);
sca10 = scales(fb10);
[~,idx60] = min(abs(sca60-1));
[~,idx10] = min(abs(sca10-1));
m60 = psi60(idx60,:);
m10 = psi10(idx10,:);

fb60 필터 뱅크의 시간과 대역폭의 곱이 더 크므로 m60 웨이블릿이 m10 웨이블릿보다 포락선 아래에 더 많은 진동을 갖는다는 사실을 확인합니다.

subplot(2,1,1)
plot(t,abs(m60))
grid on
hold on
plot(t,real(m60))
plot(t,imag(m60))
hold off
xlim([-30 30])
legend("abs(m60)","real(m60)","imag(m60)")
title("TimeBandwidth = 60")
subplot(2,1,2)
plot(t,abs(m10))
grid on
hold on
plot(t,real(m10))
plot(t,imag(m10))
hold off
xlim([-30 30])
legend("abs(m10)","real(m10)","imag(m10)")
title("TimeBandwidth = 10")

m60 크기 주파수 응답과 m10 크기 주파수 응답의 피크를 정렬합니다. m60 웨이블릿의 주파수 응답이 m10 웨이블릿의 주파수 응답보다 좁은지 확인합니다.

cf60 = centerFrequencies(fb60);
cf10 = centerFrequencies(fb10);

m60cFreq = cf60(idx60);
m10cFreq = cf10(idx10);

freqShift = 2*pi*(m60cFreq-m10cFreq);
x10 = m10.*exp(1j*freqShift*(-sigLen/2:sigLen/2-1));

figure
plot([abs(fft(m60)).' abs(fft(x10)).'])
grid on
legend("Time-bandwidth = 60","Time-bandwidth = 10")
title("Magnitude Frequency Responses")

이 예제에서는 Figure 서브플롯에 CWT 스케일로그램을 플로팅하는 방법을 보여줍니다.

음성 샘플을 불러옵니다. 데이터는 7418Hz로 샘플링됩니다. 디폴트 CWT 스케일로그램을 플로팅합니다.

load mtlb
cwt(mtlb,Fs)

신호의 연속 웨이블릿 변환과 CWT의 주파수를 구합니다.

[cfs,frq] = cwt(mtlb,Fs);

cwt 함수는 스케일로그램에 시간 좌표축과 주파수 좌표축을 설정합니다. 샘플 시간을 나타내는 벡터를 생성합니다.

tms = (0:numel(mtlb)-1)/Fs;

새 Figure에서 원래 신호를 상부 서브플롯에, 스케일로그램을 하부 서브플롯에 플로팅합니다. 주파수를 로그 스케일로 플로팅합니다.

figure
subplot(2,1,1)
plot(tms,mtlb)
axis tight
title("Signal and Scalogram")
xlabel("Time (s)")
ylabel("Amplitude")
subplot(2,1,2)
surface(tms,frq,abs(cfs))
axis tight
shading flat
xlabel("Time (s)")
ylabel("Frequency (Hz)")
set(gca,"yscale","log")