N=128;
%生成四分量信號
t1=N/4; 
t2=3*N/4; 
t3=t1; 
t4=t2; 
f1=0.15; 
f2=f1; 
f3=0.35; 
f4=f3; 
t=1:N; 
T=sqrt(2*N);
sig1=amgauss(N,t1,T).*fmconst(N,f1,t1);  
sig2=amgauss(N,t2,T).*fmconst(N,f2,t2);  
sig3=amgauss(N,t3,T).*fmconst(N,f3,t3);  
sig4=amgauss(N,t4,T).*fmconst(N,f4,t4);  
sig=sig1+sig2+sig3+sig4;
%時域波形
figure(1)
plot(1:N,real(sig));
xlabel('時間 t');
ylabel('實部幅值 A');
%頻率上下限
fmin=0.05; 
fmax=0.45; 

%計算尺度圖
[tfr,t,f]=tfrscalo(sig,t,16,fmin,fmax,N); 
Max=max(max(tfr));
V=[0.1 0.3 0.5 0.7 0.9]*Max;
figure(2)
contour(t,f,tfr,V);
xlabel('時間 t');
ylabel('頻率 f'); 
%計算affine類平滑偽Wigner分布
%%設置參數
Np=8;
Ng0=round(fliplr(linspace(2,N/10,Np/2-1)));
Nh0=round(linspace(N/4,3*N/4,Np/2-1));
[tfr,t,f]=tfrspaw(sig,t,0,Nh0(1),Ng0(1),fmin,fmax,N,1);
Max=max(max(tfr));
V=[0.1 0.3 0.5 0.7 0.9]*Max;
figure(3)
contour(t,f,tfr,V);
xlabel('時間 t');
ylabel('頻率 f'); 
%計算Wigner-Ville分布
[tfr,t,f]=tfrwv(sig);
Max=max(max(tfr)); 
V=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]*Max;
figure(4)
contour(t,f,tfr,V);
axis([1 N fmin fmax]);
xlabel('時間 t');
ylabel('頻率 f');