%樣本點數
N=4096;
%正常時心音信號
[x,fxs,bits] = wavread('normal',N);
%肺動脈瓣狹窄時心音信號
[y,fys,bits] = wavread('pulmonarystenosis',N);

%設置參數
%時域窗函數
g=window(9,'Kaiser');
%頻域窗函數
h=window(27,'Kaiser');
t=1:N;

%計算Choi-Willams分布
%%正常狀態
x=hilbert(x);
x=x-mean(x);
[tfr,t,f]=tfrcw(x,t,128,g,h,10);
%%顯示:三維圖
figure(1)
mesh(t,f*fxs/2,abs(tfr));
xlabel('樣本序號 n/個');
ylabel('頻率 f/Hz');
zlabel('幅值 A/V');
%%等高線圖
figure(2)
contour(t,f*fxs/2,abs(tfr));
xlabel('樣本序號 n/個');
ylabel('頻率 f/Hz');

%%異常狀態
y=hilbert(y);
y=y-mean(y);
[tfr,t,f]=tfrcw(y,t,128,g,h,10);

%%顯示:三維圖
figure(3)
mesh(t,f*fys/2,abs(tfr));
xlabel('樣本序號 n/個');
ylabel('頻率 f/Hz');
zlabel('幅值 A/V');

%%等高線圖
figure(4)
contour(t,f*fys/2,abs(tfr));
xlabel('樣本序號 n/個');
ylabel('頻率 f/Hz');