请问大侠们：morlet小波变换的模平方和实部的时频分布图是怎么画出来的呀？

wangyan

morlet小波变换的模平方和实部的时频分布图是怎么画出来的呀？我买的书上并没有讲这一块，请大侠们帮忙，有类似程序更是感激不尽，谢谢，谢谢，小女子将非常感谢

jasonwu

困难解决没，可以联系
QQ86914202

zhangnan3509

有什么问题可以在帖子上说呀！我们大家还能看看呢:handshake

simon21

1. function tfrview(tfr,sig,t,method,param,p1,p2,p3,p4,p5);
   
2. %TFRVIEW Visualization of time-frequency representations.
   
3. %        TFRVIEW(TFR,SIG,T,METHOD,PARAM,P1,P2,P3,P4,P5)
   
4. %        allows to visualize a time-frequency representation.
   
5. %        TFRVIEW is called through TFRQVIEW from any TFR* function.
   
6. %
   
7. %        TFR    : time-frequency representation.
   
8. %        SIG    : signal in the time-domain.  
   
9. %        T      : time instants.
   
10. %        METHOD : chosen representation (name of the corresponding M-file)
    
11. %        PARAM  : visualization parameter vector :
    
12. %         PARAM  = [DISPLAY LINLOG THRESHOLD LEVNUMB NF2 LAYOUT FS ISGRID] where
    
13. %         DISPLAY=1..5 for contour, imagesc, pcolor, surf or mesh
    
14. %         LINLOG =0/1 for linearly/logarithmically spaced levels
    
15. %                  
    
16. %         THRESHOLD  is the visualization threshold, in %
    
17. %         LEVELNUMB  is the number of levels used with contour
    
18. %         NF2        is the number of frequency bins displayed
    
19. %         LAYOUT determines the layout of the figure : TFR alone (1),  
    
20. %                TFR and SIG (2), TFR and spectrum (3), TFR and SIG and
    
21. %                spectrum (4), add 4 if you want a colorbar
    
22. %         FS     is the sampling frequency (may be set to 1.0)
    
23. %         ISGRID depends on the grids' presence :
    
24. %                isgrid=isgridsig+2*isgridspec+4*isgridtfr
    
25. %                where isgridsig=1 if a grid is present on the signal
    
26. %                and =0 if not, and so on
    
27. %         fmin   smallest normalized frequency
    
28. %         fmax   highest  normalized frequency
    
29. %        P1..P5: parameters of the representation. Run the file
    
30. %                TFRPARAM(METHOD) to know the meaning of P1..P5.           
    
31. %
    
32. %    TFRVIEW is called through TFRQVIEW by any file TFR*.
    
33. %
    
34. %            Use TFRQVIEW preferably.
    
35. %            ------------------------
    
36. %
    
37. %    See also TFRQVIEW, TFRPARAM.
    
38. 
    
39. %    F. Auger, July 1994, July 1995 -
    
40. %    O. Lemoine, October-November 1995, May-June 1996.
    
41. %       F. Auger, May 1998.
    
42. %    Copyright (c) CNRS - France 1996.
    
43. %
    
44. %    ------------------- CONFIDENTIAL PROGRAM --------------------
    
45. %    This program can not be used without the authorization of its
    
46. %    author(s). For any comment or bug report, please send e-mail to
    
47. %    f.auger@ieee.org
    
48. 
    
49. comp=computer; % so as to know the running computer
    
50. MatlabVersion=version; MatlabVersion=5;
    
51. % I hope that future Matlab versions will be more compatible
    
52. if (MatlabVersion==4),
    
53. TickLabelStr='TickLabels';
    
54. else (MatlabVersion==5),
    
55. TickLabelStr='Ticklabel';
    
56. 
    
57. end;
    
58. 
    
59. if ( nargin < 5 ),
    
60. error ('at least 5 parameters are required');
    
61. end;
    
62. 
    
63. [tfrrow,tfrcol] = size(tfr); % the size of tfr
    
64. [trow,tcol]     = size(t);   % the size of t
    
65. [Nsig,Ncol]     = size(sig); % the size of sig
    
66. maxi=max(max(tfr));
    
67. 
    
68. % Extraction of the elements of param
    
69. display     = param( 1);      % contour, imagesc, pcolor, surf, or mesh
    
70. linlog      = param( 2);      % linear or logarithmic scale
    
71. threshold   = param( 3);      % visualization threshold
    
72. levelnumb   = param( 4);      % number of levels
    
73. Nf2         = param( 5);      % number of frequency points
    
74. layout      = param( 6);      % figure layout
    
75. fs          = param( 7);      % sampling frequency
    
76. isgrid      = param( 8);      % grid(s)
    
77. fmin        = param( 9);      % smallest displayed frequency
    
78. fmax        = param(10);      % highest displayed frequency
    
79. 
    
80. if (fmin>=fmax),
    
81. error('fmin should be lower than fmax');
    
82. elseif (fmax>0.5),
    
83. error('fmax is a normalized frequency, and should be lower than 0.5');
    
84. end
    
85. 
    
86. if fs<1000, unitHz=1;                            % display in  s and  Hz
    
87. elseif (fs>=1e3 & fs<1e6), fs=fs/1e3; unitHz=2;  % display in ms and kHz
    
88. elseif (fs>=1e6), fs=fs/1e6; unitHz=3;           % display in us and MHz
    
89. end
    
90. 
    
91. linlogtfr=rem(linlog,2);
    
92. linlogspec=rem(linlog-linlogtfr,4)/2;
    
93. sigenveloppe=rem(linlog-linlogtfr-linlogspec*2,8)/4;
    
94. 
    
95. issig=rem(layout-1,2);
    
96. isspec=rem(layout-1-issig,4)/2;
    
97. iscolorbar=rem(layout-1-issig-isspec*2,8)/4;
    
98. if isempty(sig),                        % I can't make miracles
    
99. issig=0; isspec=0; layout=issig+isspec*2+1;
    
100. else  
     
101. layout=layout-4*iscolorbar;
     
102. end;
     
103. 
     
104. isgridsig =rem(isgrid,2);
     
105. isgridspec=rem(isgrid-isgridsig,4)/2;
     
106. isgridtfr =rem(isgrid-isgridsig-isgridspec*2,8)/4;
     
107. 
     
108. % Computation of isaffine and freq (vector of frequency samples)
     
109. if istfraff(method),
     
110. freq=eval(['p',num2str(nargin-5)]);       % last input argument is freqs.
     
111. isaffine=1;
     
112. if display==2,                            % imagesc not allowed
     
113.   display=3;                            % non linear scales for the axes.
     
114.   disp('Imagesc does not support non-linear scales for axes. We use pcolor instead');
     
115. end
     
116. else
     
117. isaffine=0;                            % weyl-heisenberg group of distributions
     
118. freq=(0.5*(0:Nf2-1)/Nf2);              % dispaly only the positive frequencies
     
119. end
     
120. freqr=freq*fs; ts=t/fs;                 % real time and frequency
     
121. 
     
122. 
     
123. % Update variables mini, levels, LinLogStr according to linlog
     
124. if ~linlogtfr,
     
125. if (display==4)|(display==5),          % surf or mesh
     
126.   mini=min(min(tfr));
     
127. else
     
128.   mini=max(min(min(tfr)),maxi*threshold/100.0);
     
129. end
     
130. levels=linspace(mini,maxi,levelnumb+1);
     
131. LinLogStr=', lin. scale';
     
132. else
     
133. mini=max(min(min(tfr)),maxi*threshold/100.0);
     
134. levels=logspace(log10(mini),log10(maxi),levelnumb+1);
     
135. LinLogStr=', log. scale';
     
136. end;
     
137. 
     
138. % Initialization of the current figure
     
139. zoom off; clf;
     
140. set(gcf,'Resize','On','NextPlot','Add');
     
141. 
     
142. % Initialization of the axes
     
143. if iscolorbar,
     
144. axcb   = axes('Units','normal','Visible','off','Box','On');
     
145. set(gcf,'UserData',[get(gcf,'UserData') axcb]);
     
146. end;
     
147. 
     
148. if issig,
     
149. axsig  = axes('Units','normal','Visible','off','Box','On');
     
150. if comp(1:2)=='PC', set(axsig ,'fontsize',10); end
     
151. set(gcf,'UserData',[get(gcf,'UserData') axsig]);
     
152. end;
     
153. 
     
154. if isspec,
     
155. axspec = axes('Units','normal','Visible','off','Box','On');
     
156. if comp(1:2)=='PC', set(axspec,'fontsize',10); end;
     
157. set(gcf,'UserData',[get(gcf,'UserData') axspec]);
     
158. end;
     
159. 
     
160. axtfr  = axes('Units','normal','Visible','off','Box','On');
     
161. if comp(1:2)=='PC', set(axtfr ,'fontsize',10); end
     
162. set(gcf,'UserData',[get(gcf,'UserData') axtfr]);
     
163. 
     
164. % Test of analycity and computation of spec
     
165. if ~isempty(sig),
     
166.   for k=1:Ncol,
     
167.    isana=1; alpha=2; Lt=max(t)-min(t)+1;      
     
168.    if 2*Nf2>=Lt,
     
169.     spec(:,k)=abs(fft(sig(min(t):max(t),k),2*Nf2)).^2;
     
170.    else
     
171.     % modifications :  F. Auger (fog), 30/11/97
     
172.     Nb_tranches_fog = floor(Lt/(2*Nf2));
     
173.     % fprintf('%f \n',Nb_tranches_fog);
     
174.     spec(:,k)=zeros(2*Nf2,1);
     
175.     for Num_tranche_fog=0:Nb_tranches_fog-1,
     
176.      % fprintf('%f \n',Num_tranche_fog);
     
177.      spec(:,k)=spec(:,k)+abs(fft(sig(min(t)+2*Nf2*Num_tranche_fog+(0:2*Nf2-1),k))).^2;
     
178.     end;
     
179. 
     
180.     if (Lt>Nb_tranches_fog*2*Nf2),
     
181.      spectre_fog=fft(sig(min(t)+tfrrow*Nb_tranches_fog:max(t),k),2*Nf2);
     
182.      spectre_fog=spectre_fog(:);
     
183.      spec(:,k)=spec(:,k)+abs(spectre_fog).^2;
     
184.     end;   
     
185.     % sp=abs(fft(sig(min(t):max(t),k))).^2;
     
186.     % fr1=(0.5*(0:Lt-1)/Lt)*fs;
     
187.     % fr2=(0.5*(0:2*Nf2-1)/2/Nf2)*fs;
     
188.     % spec(:,k)=interp1(fr1,sp,fr2);
     
189.    end
     
190.    spec1=sum(spec(1:Nf2,k));
     
191.    spec2=sum(spec(Nf2+1:2*Nf2,k));
     
192.    if spec2>0.1*spec1,
     
193.     isana=0;
     
194.     if ~isreal(sig(min(t):max(t),k)),
     
195.      alpha=1;
     
196.     end
     
197.    end
     
198.   end
     
199. end
     
200. 
     
201. if layout==1,                          % Time-Frequency Representation only
     
202.   set(axtfr,'Position',[0.10 0.10 0.80 0.80]);
     
203. 
     
204. elseif layout==2,            % TFR + Signal
     
205.   set(axtfr,'Position',[0.10 0.10 0.80 0.55]);
     
206.   set(axsig,'Position',[0.10 0.73 0.80 0.20]);
     
207.   axes(axsig);  
     
208.   if sigenveloppe,
     
209.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)),...
     
210.         (min(t):max(t))/fs, abs(sig(min(t):max(t),:)));
     
211.   else
     
212.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)));
     
213.   end;
     
214.   set(gca,['X' TickLabelStr],[]);
     
215.   ylabel('Real part');
     
216.   title('Signal in time');
     
217.   Min=min(min(real(sig))); Max=max(max(real(sig)));
     
218.   axis([min(ts) max(ts) Min Max]);
     
219. 
     
220. elseif layout==3,            % TFR + spectrum
     
221.   set(axspec,'Position',[0.10 0.10 0.15 0.80]);
     
222.   set(axtfr ,'Position',[0.35 0.10 0.55 0.80]);
     
223. 
     
224.   axes(axspec);
     
225.   if isaffine,
     
226.    f1=freqr(1); f2=freqr(Nf2); df=f2-f1;
     
227.    Nf4=round((Nf2-1)*fs/(2*df))+1;
     
228.    for k=1:Ncol,
     
229.     spec(1:alpha*Nf4,k)=abs(fft(sig(min(t):max(t),k),alpha*Nf4)).^2;
     
230.    end
     
231.    spec=spec((round(f1*2*(Nf4-1)/fs)+1):(round(f1*2*(Nf4-1)/fs)+Nf2),:);
     
232.    freqs=linspace(f1,f2,Nf2);
     
233.   else
     
234.    freqs=freqr;
     
235.    spec=spec(1:Nf2,:);
     
236.   end
     
237.   Maxsp=max(max(spec));
     
238.   if linlogspec==0,
     
239.    plot(freqs,spec);
     
240.    title('Linear scale');
     
241.    % set(axspec,'ytick',[0 Maxsp*max(eps,threshold)*0.01 Maxsp]);
     
242.    set(axspec,'YTickMode', 'auto');
     
243.    set(axspec,'Ylim', [Maxsp*threshold*0.01 Maxsp*1.2]);
     
244.    set(axspec,'Xlim', [fmin fmax]);
     
245.   else
     
246.    plot(freqs,10*log10(spec/Maxsp));
     
247.    title('Log. scale [dB]');
     
248.    set(axspec,'YTickMode', 'auto');
     
249.    set(axspec,'Ylim',[10*log10(threshold*0.01) 0]);
     
250.    set(axspec,'Xlim', [fmin fmax]);
     
251.   end
     
252.   xlabel('Energy spectral density');
     
253.   Nsp=length(spec);
     
254.   set(gca, ['X' TickLabelStr],[],'view',[-90 90]);
     
255. 
     
256. elseif layout==4,            % TFR + signal + spectrum
     
257. 
     
258.   set(axspec,'Position',[0.10 0.10 0.15 0.55]);
     
259.   set(axsig ,'Position',[0.35 0.73 0.55 0.20]);
     
260.   set(axtfr ,'Position',[0.35 0.10 0.55 0.55]);
     
261. 
     
262.   axes(axsig);
     
263.   if sigenveloppe,
     
264.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)),...
     
265.         (min(t):max(t))/fs, abs(sig(min(t):max(t),:)));
     
266.   else
     
267.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)));
     
268.   end;
     
269.   ylabel('Real part');
     
270.   title('Signal in time');
     
271.   set(gca,['X' TickLabelStr],[]);
     
272. 
     
273.   Min=min(min(real(sig))); Max=max(max(real(sig)));
     
274.   axis([min(ts) max(ts) Min Max]);
     
275. 
     
276.   axes(axspec);
     
277.   if isaffine,
     
278. % IF YOU UNDERSTAND SOMETHING TO THAT, PLEASE EXPLAIN ME (f. auger)
     
279. %   f1=freqr(1); f2=freqr(Nf2); df=f2-f1;
     
280. %   Nf4=round((Nf2-1)*fs/(2*df))+1;
     
281. %   for k=1:Ncol,
     
282. %    spec(1:alpha*Nf4,k)=abs(fft(sig(min(t):max(t),k),alpha*Nf4)).^2;
     
283. %   end
     
284. %   spec=spec((round(f1*2*(Nf4-1)/fs)+1):(round(f1*2*(Nf4-1)/fs)+Nf2),:);
     
285. %   freqs=linspace(f1,f2,Nf2);
     
286.    for k=1:Ncol,
     
287.     freqs=linspace(freqr(1),freqr(Nf2),Nf2);
     
288.     spec=interp1(0.5*fs*(0:Nf2-1)/Nf2,spec(1:Nf2,k),freqs);
     
289.    end;
     
290.   else
     
291.    freqs=freqr;
     
292.    spec=spec(1:Nf2,:);
     
293.   end
     
294.   Maxsp=max(max(spec));
     
295.   if linlogspec==0,
     
296.    plot(freqs,spec);
     
297.    title('Linear scale');  
     
298.    set(axspec,'YTickMode', 'auto');
     
299.    set(axspec,'Ylim', [Maxsp*threshold*0.01 Maxsp*1.2]);
     
300.    set(axspec,'Xlim', [fmin*fs fmax*fs]);
     
301.   else
     
302.    plot(freqs,10*log10(spec/Maxsp));
     
303.    title('Log. scale [dB]');
     
304.    set(axspec,'Ytickmode','auto');
     
305.    set(axspec,'Ylim',[10*log10(threshold*0.01) 0]);
     
306.    set(axspec,'Xlim', [fmin*fs fmax*fs]);
     
307.   end
     
308.   xlabel('Energy spectral density');
     
309.   Nsp=length(spec);
     
310.   set(gca,['X' TickLabelStr],[],'view',[-90 90]);
     
311. 
     
312. end;
     
313. 
     
314.   if iscolorbar,                     % Is there a color bar ?
     
315.    PositionTfr=get(axtfr,'Position');
     
316.    set(axtfr,'Position',PositionTfr-[0 0 0.03 0]);
     
317.    set(axcb, 'Position',[PositionTfr(1)+PositionTfr(3)-0.01,...
     
318.                          PositionTfr(2) 0.01 PositionTfr(4)]);
     
319.    axes(axcb);
     
320.    Ncolors=length(colormap);
     
321.    [cmin,cmax]=caxis; colorvect=linspace(cmax,cmin,Ncolors);
     
322.    imagesc(colorvect'); axis('off');
     
323.    if issig,                          % there is a signal
     
324.     PositionSig=get(axsig,'Position');
     
325.     set(axsig,'Position',PositionSig-[0 0 0.03 0]);
     
326.    end
     
327.   end
     
328. 
     
329. axes(axtfr);                         % Display the tfr
     
330. if (tcol==1),
     
331.   plot(freqr,tfr(1:Nf2));
     
332.   set(axtfr,'Xlim', [fmin*fs fmax*fs]);
     
333. 
     
334. else
     
335.   if strcmp(computer,'MAC'),
     
336.    tfr=flipud(tfr(1:Nf2,:));
     
337.   else
     
338.    tfr=tfr(1:Nf2,:);
     
339.   end;
     
340. 
     
341.   indmin=find(tfr<mini);
     
342.   tfr(indmin)=mini*ones(1,length(indmin));
     
343. 
     
344.   indmax=find(tfr>maxi);
     
345.   tfr(indmax)=maxi*ones(1,length(indmax));
     
346.   if isaffine & (display==2),
     
347.    fprintf('imagesc does not support non-linear scales for axes. Replaced by pcolor.\n');
     
348.    display=3;
     
349.   end;
     
350. 
     
351.   if display==1,                  % contour
     
352.    contour(ts,freqr,tfr,levels);  % contour(tfr,levels,ts,freqr);
     
353.    set(axtfr,'Ylim', [fmin*fs fmax*fs]);
     
354.    DisplayStr=', contour';
     
355.   elseif display==2,              % imagesc
     
356.    if linlogtfr==0,
     
357.     imagesc(ts,freqr,tfr); axis('xy');
     
358.    else
     
359.     imagesc(ts,freqr,log10(tfr));axis('xy');
     
360.    end
     
361.    set(axtfr,'Ylim', [fmin*fs fmax*fs]);
     
362.    DisplayStr=', imagesc';
     
363. 
     
364.   elseif display==3,              % pcolor
     
365.    if linlogtfr==0,
     
366.     pcolor(ts,freqr,tfr); shading interp;
     
367.    else
     
368.     pcolor(ts,freqr,log10(tfr)); shading interp;
     
369.    end
     
370.    set(axtfr,'Ylim', [fmin*fs fmax*fs]);
     
371.    DisplayStr=', pcolor';
     
372. 
     
373.   elseif display==4,              % surf
     
374.    if linlogtfr==0,
     
375.     surf(ts,freqr,tfr); shading interp;
     
376.     zlabel('Amplitude');
     
377.     axis([ts(1) ts(tcol) fmin*fs fmax*fs mini maxi]);
     
378.    else
     
379.     surf(ts,freqr,log10(tfr)); shading interp;
     
380.     zlabel('Positive values');
     
381.     axis([ts(1) ts(tcol) fmin fmax log10(mini) log10(maxi)]);
     
382.    end
     
383.    DisplayStr=', surf';
     
384. 
     
385.   elseif display==5,              % mesh
     
386.    if linlogtfr==0,
     
387.     mesh(ts,freqr,tfr); shading interp;
     
388.     zlabel('Amplitude');
     
389.     axis([ts(1) ts(tcol) fmin*fs fmax*fs mini maxi]);
     
390.    else
     
391.     mesh(ts,freqr,log10(tfr));shading interp;
     
392.     zlabel('Positive values');
     
393.     axis([ts(1) ts(tcol) fmin*fs fmax*fs log10(mini) log10(maxi)]);
     
394.    end
     
395.    DisplayStr=', mesh';
     
396.   end
     
397. 
     
398. % Define the title and check the input arguments depending on 'method'
     
399. 
     
400. method=method(4:length(method));
     
401. 
     
402. if nargin==5, % if there is no additional parameters, do the best.
     
403.   title([method, LinLogStr,DisplayStr,...
     
404.          ', Threshold=',num2str(threshold),'%']);
     
405. 
     
406. elseif strcmp(method,'WV'  ) | strcmp(method,'MH') | ...
     
407.         strcmp(method,'PAGE'), % no parameters
     
408.   title([method,', Nf=',num2str(Nf2), LinLogStr,DisplayStr,...
     
409.         ', Threshold=',num2str(threshold),'%']);
     
410. 
     
411. elseif strcmp(method,'PWV'  )|strcmp(method,'PMH'  )| ...
     
412.         strcmp(method,'SP'   )|strcmp(method,'PPAGE')| ...
     
413.         strcmp(method,'RSP'  )|strcmp(method,'RPPAG')| ...
     
414.         strcmp(method,'RPWV' )|strcmp(method,'RPMH' ),
     
415.   h=p1;[hrow,hcol]=size(h); Lh=(hrow-1)/2; % one parameter
     
416.   if (hcol~=1)|(rem(hrow,2)==0),
     
417.    error('h must be a smoothing window with odd length'); end;
     
418.   title([method, ', Lh=',num2str(Lh), ', Nf=',num2str(Nf2),...
     
419.         LinLogStr, DisplayStr,', Threshold=',num2str(threshold),'%']);
     
420. 
     
421. elseif strcmp(method,'STFT'), % short-time fourier transform case
     
422.   h=p1;[hrow,hcol]=size(h); Lh=(hrow-1)/2;
     
423.   if (hcol~=1)|(rem(hrow,2)==0),
     
424.    error('h must be a smoothing window with odd length'); end;
     
425.   title(['|',method, '|^2, Lh=',num2str(Lh),...
     
426.         ', Nf=',num2str(Nf2), LinLogStr, DisplayStr,', Thld=',...
     
427.         num2str(threshold),'%']);
     
428. 
     
429. elseif strcmp(method,'SPWV' ) | strcmp(method,'MHS'  )| ...
     
430.         strcmp(method,'RSPWV') | strcmp(method,'RMHS' )| ...
     
431.         strcmp(method,'ZAM'  ) | strcmp(method,'RIDBN')|...
     
432.         strcmp(method,'BJ'   ) | strcmp(method,'RIDB' )| ...
     
433.         strcmp(method,'RIDH' ) | strcmp(method,'RIDT' ),
     
434.   g=p1; [grow,gcol]=size(g); Lg=(grow-1)/2;
     
435.   if (gcol~=1)|(rem(grow,2)==0),
     
436.    error('g must be a smoothing window with odd length'); end;
     
437.   h=p2; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
     
438.   if (hcol~=1)|(rem(hrow,2)==0),
     
439.    error('h must be a smoothing window with odd length'); end;
     
440.   title([method,', Lg=',num2str(Lg),', Lh=',num2str(Lh),...
     
441.          ', Nf=',num2str(Nf2),LinLogStr, DisplayStr,...
     
442.          ', Threshold=',num2str(threshold),'%']);
     
443. 
     
444. elseif strcmp(method,'MMCE'),
     
445.   h=p1;[hrow,hcol]=size(h); Lh=(hrow-1)/2;
     
446.   if (rem(hrow,2)==0),
     
447.    error('h must be a smoothing window with odd length'); end;
     
448.   title([method, ', Lh=',num2str(Lh), ', Nf=',num2str(Nf2),...
     
449.         LinLogStr, DisplayStr,', Threshold=',num2str(threshold),'%']);
     
450. 
     
451. elseif strcmp(method,'CW' ) | strcmp(method,'BUD'),
     
452.   g=p1; [grow,gcol]=size(g); Lg=(grow-1)/2;
     
453.   if (gcol~=1)|(rem(grow,2)==0),
     
454.    error('g must be a smoothing window with odd length'); end;
     
455.   h=p2; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
     
456.   if (hcol~=1)|(rem(hrow,2)==0),
     
457.    error('h must be a smoothing window with odd length'); end;
     
458.   sigma=p3;
     
459.   title([method,', Lg=',num2str(Lg),', Lh=',num2str(Lh),...
     
460.          ' sigma=',num2str(sigma),', Nf=',num2str(Nf2),...
     
461.          LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);
     
462. 
     
463. elseif strcmp(method,'GRD')
     
464.   g=p1; [grow,gcol]=size(g); Lg=(grow-1)/2;
     
465.   if (gcol~=1)|(rem(grow,2)==0),
     
466.    error('g must be a smoothing window with odd length'); end;
     
467.   h=p2; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
     
468.   if (hcol~=1)|(rem(hrow,2)==0),
     
469.    error('h must be a smoothing window with odd length'); end;
     
470.   title([method,', Lg=',num2str(Lg),', Lh=',num2str(Lh),...
     
471.          ', rs =',num2str(p3), ', M/N =',num2str(p4), ...
     
472.          ', Nf =',num2str(Nf2), ...
     
473.          LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);
     
474. 
     
475. elseif strcmp(method,'MSC' ) | strcmp(method,'RMSC' )
     
476.   f0T=p1; if (f0T<=0), error('f0T must be positive'); end;
     
477.   title([method, ', f0T=',num2str(f0T), ', Nf=',num2str(Nf2),...
     
478.         LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);
     
479. 
     
480. elseif strcmp(method,'RGAB' )
     
481.   Nh=p1; if (Nh<=0), error('Nh must be positive'); end;
     
482.   title([method, ', Nh=',num2str(Nh), ', Nf=',num2str(Nf2),...
     
483.         LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);
     
484. 
     
485. 
     
486. elseif strcmp(method,'DFLA' ) | strcmp(method,'UNTER' )| ...
     
487.         strcmp(method,'BERT' ),
     
488.   N=p1;
     
489.   if (N<=0),                error('N must be positive'); end;
     
490.   title([method, ', N=',num2str(N), LinLogStr, DisplayStr, ', Threshold=',...
     
491.          num2str(threshold), '%']);  
     
492. 
     
493. elseif strcmp(method,'SCALO'),
     
494.   Nh0=p1; N=p2;
     
495.   if (Nh0<0),               error('Nh0 must be positive'); end;
     
496.   if (N<=0),                error('N must be positive'); end;
     
497.   if (Nh0>0),
     
498.    title([method, ', Morlet wavelet, Nh0=', num2str(Nh0), ...
     
499.          ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
     
500.          num2str(threshold), '%']);  
     
501.   else
     
502.    title([method, ', Mexican hat, N=',num2str(N), LinLogStr, DisplayStr, ...
     
503.           ', Thld=', num2str(threshold), '%']);  
     
504.   end
     
505. 
     
506. elseif strcmp(method,'ASPW'),
     
507.   Nh0=p1; Ng0=p2; N=p3;
     
508.   if (Nh0<0),               error('Nh0 must be positive'); end;
     
509.   if (Ng0<0),               error('Ng0 must be positive'); end;
     
510.   if (N<=0),                error('N must be positive'); end;
     
511.   if (Nh0>0),
     
512.    title([method, ', Morlet wlt, Nh0=', num2str(Nh0), ', Ng0=',...
     
513.          num2str(Ng0), ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
     
514.          num2str(threshold), '%']);  
     
515.   else
     
516.    title([method, ', Mexican hat, Ng0=',num2str(Ng0),...  
     
517.          ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
     
518.          num2str(threshold), '%']);  
     
519.   end
     
520. 
     
521. elseif strcmp(method,'SPAW'),
     
522.   K=p1; Nh0=p2; Ng0=p3; N=p4;
     
523.   if (Nh0<0),               error('Nh0 must be positive'); end;
     
524.   if (Ng0<0),               error('Ng0 must be positive'); end;
     
525.   if (N<=0),                error('N must be positive'); end;
     
526.   if (Nh0>0),
     
527.    title([method, ', K=', num2str(K), ', Morlet wlt, Nh0=',...
     
528.          num2str(Nh0), ', Ng0=',...
     
529.          num2str(Ng0), ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
     
530.          num2str(threshold), '%']);  
     
531.   else
     
532.    title([method, ', K=', num2str(K), ', Mexican hat, Ng0=',...
     
533.          num2str(Ng0),', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
     
534.          num2str(threshold), '%']);  
     
535.   end
     
536. 
     
537. 
     
538. elseif strcmp(method,'GABOR'),
     
539.   N=p1; Q=p2; h=p3; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
     
540.   if (hcol~=1)|(rem(hrow,2)==0),
     
541.    error('h must be a smoothing window with odd length'); end;
     
542.   title([method, ', Lh=',num2str(Lh), ', Nf=',...
     
543.          num2str(Nf2),', N=',num2str(N),', Q=',num2str(Q),...
     
544.          LinLogStr, DisplayStr, ', Thld=',num2str(threshold),'%']);
     
545. 
     
546. end;
     
547. end
     
548. 
     
549. % add the correct legend on the axes
     
550. if unitHz==1,
     
551. xlabel('Time [s]'); ylabel('Frequency [Hz]');
     
552. elseif unitHz==2,
     
553. xlabel('Time [ms]'); ylabel('Frequency [kHz]');
     
554. elseif unitHz==3,
     
555. xlabel('Time [祍]'); ylabel('Frequency [MHz]');
     
556. end
     
557. 
     
558. 
     
559. if (isgridsig & issig),        % Updating of the grids
     
560. axes(axsig); grid on   
     
561. elseif (~isgridsig & issig),
     
562. axes(axsig); grid off
     
563. end
     
564. if (isgridspec & isspec),
     
565. axes(axspec); grid on
     
566. elseif (~isgridspec & isspec),
     
567. axes(axspec); grid off
     
568. end
     
569. 
     
570. if (isgridtfr),            % upating of the grid on the tfr
     
571. axes(axtfr); grid on  
     
572. elseif (~isgridtfr),
     
573. axes(axtfr); grid off
     
574. end

复制代码

jasonwu

就是我在matlab中输入命令：coefs=cwt(r,scales,'cmor 1-1.5')，得到小波变换系数，是一个81*32的矩阵，然后求其实部和模平方分别为：x=real(coefs),y=abs(coefs).^2，同样x,y也是81*32的矩阵，然后画等值线图：contour(x),contour(y)。
    问题1就是，在matlab中画出来的等直线图，我不会调疏密，还有就是值小于0的线调成虚线，还有上式中1-1.5是什么意思啊？
    问题2就是，我今天尝试用surfer来画等直线图，但我不知道怎么把x一个81*32的矩阵转化成surfer合法的输入，还请大侠不吝赐教。

oppo

有谁知道?感谢ing~~~~~~~~~~

stao136

也有同样的疑问

stao136

多学习学习？小波方差是怎么做的？

fudy10

回复 jasonwu 的帖子
建议你在版内搜索“连续小波变换”