声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4416|回复: 8

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

[复制链接]
发表于 2006-7-24 11:13 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

使用道具 举报

发表于 2007-4-27 14:39 | 显示全部楼层

困难解决没有?

困难解决没,可以联系
QQ86914202
发表于 2007-4-27 14:46 | 显示全部楼层

回复 #2 jasonwu 的帖子

有什么问题可以在帖子上说呀!我们大家还能看看呢:handshake
发表于 2007-4-27 14:51 | 显示全部楼层
  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. %    F. Auger, July 1994, July 1995 -
  39. %    O. Lemoine, October-November 1995, May-June 1996.
  40. %       F. Auger, May 1998.
  41. %    Copyright (c) CNRS - France 1996.
  42. %
  43. %    ------------------- CONFIDENTIAL PROGRAM --------------------
  44. %    This program can not be used without the authorization of its
  45. %    author(s). For any comment or bug report, please send e-mail to
  46. %    f.auger@ieee.org

  47. comp=computer; % so as to know the running computer
  48. MatlabVersion=version; MatlabVersion=5;
  49. % I hope that future Matlab versions will be more compatible
  50. if (MatlabVersion==4),
  51. TickLabelStr='TickLabels';
  52. else (MatlabVersion==5),
  53. TickLabelStr='Ticklabel';

  54. end;

  55. if ( nargin < 5 ),
  56. error ('at least 5 parameters are required');
  57. end;

  58. [tfrrow,tfrcol] = size(tfr); % the size of tfr
  59. [trow,tcol]     = size(t);   % the size of t
  60. [Nsig,Ncol]     = size(sig); % the size of sig
  61. maxi=max(max(tfr));

  62. % Extraction of the elements of param
  63. display     = param( 1);      % contour, imagesc, pcolor, surf, or mesh
  64. linlog      = param( 2);      % linear or logarithmic scale
  65. threshold   = param( 3);      % visualization threshold
  66. levelnumb   = param( 4);      % number of levels
  67. Nf2         = param( 5);      % number of frequency points
  68. layout      = param( 6);      % figure layout
  69. fs          = param( 7);      % sampling frequency
  70. isgrid      = param( 8);      % grid(s)
  71. fmin        = param( 9);      % smallest displayed frequency
  72. fmax        = param(10);      % highest displayed frequency

  73. if (fmin>=fmax),
  74. error('fmin should be lower than fmax');
  75. elseif (fmax>0.5),
  76. error('fmax is a normalized frequency, and should be lower than 0.5');
  77. end

  78. if fs<1000, unitHz=1;                            % display in  s and  Hz
  79. elseif (fs>=1e3 & fs<1e6), fs=fs/1e3; unitHz=2;  % display in ms and kHz
  80. elseif (fs>=1e6), fs=fs/1e6; unitHz=3;           % display in us and MHz
  81. end

  82. linlogtfr=rem(linlog,2);
  83. linlogspec=rem(linlog-linlogtfr,4)/2;
  84. sigenveloppe=rem(linlog-linlogtfr-linlogspec*2,8)/4;

  85. issig=rem(layout-1,2);
  86. isspec=rem(layout-1-issig,4)/2;
  87. iscolorbar=rem(layout-1-issig-isspec*2,8)/4;
  88. if isempty(sig),                        % I can't make miracles
  89. issig=0; isspec=0; layout=issig+isspec*2+1;
  90. else  
  91. layout=layout-4*iscolorbar;
  92. end;

  93. isgridsig =rem(isgrid,2);
  94. isgridspec=rem(isgrid-isgridsig,4)/2;
  95. isgridtfr =rem(isgrid-isgridsig-isgridspec*2,8)/4;

  96. % Computation of isaffine and freq (vector of frequency samples)
  97. if istfraff(method),
  98. freq=eval(['p',num2str(nargin-5)]);       % last input argument is freqs.
  99. isaffine=1;
  100. if display==2,                            % imagesc not allowed
  101.   display=3;                            % non linear scales for the axes.
  102.   disp('Imagesc does not support non-linear scales for axes. We use pcolor instead');
  103. end
  104. else
  105. isaffine=0;                            % weyl-heisenberg group of distributions
  106. freq=(0.5*(0:Nf2-1)/Nf2);              % dispaly only the positive frequencies
  107. end
  108. freqr=freq*fs; ts=t/fs;                 % real time and frequency


  109. % Update variables mini, levels, LinLogStr according to linlog
  110. if ~linlogtfr,
  111. if (display==4)|(display==5),          % surf or mesh
  112.   mini=min(min(tfr));
  113. else
  114.   mini=max(min(min(tfr)),maxi*threshold/100.0);
  115. end
  116. levels=linspace(mini,maxi,levelnumb+1);
  117. LinLogStr=', lin. scale';
  118. else
  119. mini=max(min(min(tfr)),maxi*threshold/100.0);
  120. levels=logspace(log10(mini),log10(maxi),levelnumb+1);
  121. LinLogStr=', log. scale';
  122. end;

  123. % Initialization of the current figure
  124. zoom off; clf;
  125. set(gcf,'Resize','On','NextPlot','Add');

  126. % Initialization of the axes
  127. if iscolorbar,
  128. axcb   = axes('Units','normal','Visible','off','Box','On');
  129. set(gcf,'UserData',[get(gcf,'UserData') axcb]);
  130. end;

  131. if issig,
  132. axsig  = axes('Units','normal','Visible','off','Box','On');
  133. if comp(1:2)=='PC', set(axsig ,'fontsize',10); end
  134. set(gcf,'UserData',[get(gcf,'UserData') axsig]);
  135. end;

  136. if isspec,
  137. axspec = axes('Units','normal','Visible','off','Box','On');
  138. if comp(1:2)=='PC', set(axspec,'fontsize',10); end;
  139. set(gcf,'UserData',[get(gcf,'UserData') axspec]);
  140. end;

  141. axtfr  = axes('Units','normal','Visible','off','Box','On');
  142. if comp(1:2)=='PC', set(axtfr ,'fontsize',10); end
  143. set(gcf,'UserData',[get(gcf,'UserData') axtfr]);

  144. % Test of analycity and computation of spec
  145. if ~isempty(sig),
  146.   for k=1:Ncol,
  147.    isana=1; alpha=2; Lt=max(t)-min(t)+1;      
  148.    if 2*Nf2>=Lt,
  149.     spec(:,k)=abs(fft(sig(min(t):max(t),k),2*Nf2)).^2;
  150.    else
  151.     % modifications :  F. Auger (fog), 30/11/97
  152.     Nb_tranches_fog = floor(Lt/(2*Nf2));
  153.     % fprintf('%f \n',Nb_tranches_fog);
  154.     spec(:,k)=zeros(2*Nf2,1);
  155.     for Num_tranche_fog=0:Nb_tranches_fog-1,
  156.      % fprintf('%f \n',Num_tranche_fog);
  157.      spec(:,k)=spec(:,k)+abs(fft(sig(min(t)+2*Nf2*Num_tranche_fog+(0:2*Nf2-1),k))).^2;
  158.     end;

  159.     if (Lt>Nb_tranches_fog*2*Nf2),
  160.      spectre_fog=fft(sig(min(t)+tfrrow*Nb_tranches_fog:max(t),k),2*Nf2);
  161.      spectre_fog=spectre_fog(:);
  162.      spec(:,k)=spec(:,k)+abs(spectre_fog).^2;
  163.     end;   
  164.     % sp=abs(fft(sig(min(t):max(t),k))).^2;
  165.     % fr1=(0.5*(0:Lt-1)/Lt)*fs;
  166.     % fr2=(0.5*(0:2*Nf2-1)/2/Nf2)*fs;
  167.     % spec(:,k)=interp1(fr1,sp,fr2);
  168.    end
  169.    spec1=sum(spec(1:Nf2,k));
  170.    spec2=sum(spec(Nf2+1:2*Nf2,k));
  171.    if spec2>0.1*spec1,
  172.     isana=0;
  173.     if ~isreal(sig(min(t):max(t),k)),
  174.      alpha=1;
  175.     end
  176.    end
  177.   end
  178. end

  179. if layout==1,                          % Time-Frequency Representation only
  180.   set(axtfr,'Position',[0.10 0.10 0.80 0.80]);

  181. elseif layout==2,            % TFR + Signal
  182.   set(axtfr,'Position',[0.10 0.10 0.80 0.55]);
  183.   set(axsig,'Position',[0.10 0.73 0.80 0.20]);
  184.   axes(axsig);  
  185.   if sigenveloppe,
  186.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)),...
  187.         (min(t):max(t))/fs, abs(sig(min(t):max(t),:)));
  188.   else
  189.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)));
  190.   end;
  191.   set(gca,['X' TickLabelStr],[]);
  192.   ylabel('Real part');
  193.   title('Signal in time');
  194.   Min=min(min(real(sig))); Max=max(max(real(sig)));
  195.   axis([min(ts) max(ts) Min Max]);

  196. elseif layout==3,            % TFR + spectrum
  197.   set(axspec,'Position',[0.10 0.10 0.15 0.80]);
  198.   set(axtfr ,'Position',[0.35 0.10 0.55 0.80]);

  199.   axes(axspec);
  200.   if isaffine,
  201.    f1=freqr(1); f2=freqr(Nf2); df=f2-f1;
  202.    Nf4=round((Nf2-1)*fs/(2*df))+1;
  203.    for k=1:Ncol,
  204.     spec(1:alpha*Nf4,k)=abs(fft(sig(min(t):max(t),k),alpha*Nf4)).^2;
  205.    end
  206.    spec=spec((round(f1*2*(Nf4-1)/fs)+1):(round(f1*2*(Nf4-1)/fs)+Nf2),:);
  207.    freqs=linspace(f1,f2,Nf2);
  208.   else
  209.    freqs=freqr;
  210.    spec=spec(1:Nf2,:);
  211.   end
  212.   Maxsp=max(max(spec));
  213.   if linlogspec==0,
  214.    plot(freqs,spec);
  215.    title('Linear scale');
  216.    % set(axspec,'ytick',[0 Maxsp*max(eps,threshold)*0.01 Maxsp]);
  217.    set(axspec,'YTickMode', 'auto');
  218.    set(axspec,'Ylim', [Maxsp*threshold*0.01 Maxsp*1.2]);
  219.    set(axspec,'Xlim', [fmin fmax]);
  220.   else
  221.    plot(freqs,10*log10(spec/Maxsp));
  222.    title('Log. scale [dB]');
  223.    set(axspec,'YTickMode', 'auto');
  224.    set(axspec,'Ylim',[10*log10(threshold*0.01) 0]);
  225.    set(axspec,'Xlim', [fmin fmax]);
  226.   end
  227.   xlabel('Energy spectral density');
  228.   Nsp=length(spec);
  229.   set(gca, ['X' TickLabelStr],[],'view',[-90 90]);

  230. elseif layout==4,            % TFR + signal + spectrum

  231.   set(axspec,'Position',[0.10 0.10 0.15 0.55]);
  232.   set(axsig ,'Position',[0.35 0.73 0.55 0.20]);
  233.   set(axtfr ,'Position',[0.35 0.10 0.55 0.55]);

  234.   axes(axsig);
  235.   if sigenveloppe,
  236.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)),...
  237.         (min(t):max(t))/fs, abs(sig(min(t):max(t),:)));
  238.   else
  239.    plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)));
  240.   end;
  241.   ylabel('Real part');
  242.   title('Signal in time');
  243.   set(gca,['X' TickLabelStr],[]);

  244.   Min=min(min(real(sig))); Max=max(max(real(sig)));
  245.   axis([min(ts) max(ts) Min Max]);

  246.   axes(axspec);
  247.   if isaffine,
  248. % IF YOU UNDERSTAND SOMETHING TO THAT, PLEASE EXPLAIN ME (f. auger)
  249. %   f1=freqr(1); f2=freqr(Nf2); df=f2-f1;
  250. %   Nf4=round((Nf2-1)*fs/(2*df))+1;
  251. %   for k=1:Ncol,
  252. %    spec(1:alpha*Nf4,k)=abs(fft(sig(min(t):max(t),k),alpha*Nf4)).^2;
  253. %   end
  254. %   spec=spec((round(f1*2*(Nf4-1)/fs)+1):(round(f1*2*(Nf4-1)/fs)+Nf2),:);
  255. %   freqs=linspace(f1,f2,Nf2);
  256.    for k=1:Ncol,
  257.     freqs=linspace(freqr(1),freqr(Nf2),Nf2);
  258.     spec=interp1(0.5*fs*(0:Nf2-1)/Nf2,spec(1:Nf2,k),freqs);
  259.    end;
  260.   else
  261.    freqs=freqr;
  262.    spec=spec(1:Nf2,:);
  263.   end
  264.   Maxsp=max(max(spec));
  265.   if linlogspec==0,
  266.    plot(freqs,spec);
  267.    title('Linear scale');  
  268.    set(axspec,'YTickMode', 'auto');
  269.    set(axspec,'Ylim', [Maxsp*threshold*0.01 Maxsp*1.2]);
  270.    set(axspec,'Xlim', [fmin*fs fmax*fs]);
  271.   else
  272.    plot(freqs,10*log10(spec/Maxsp));
  273.    title('Log. scale [dB]');
  274.    set(axspec,'Ytickmode','auto');
  275.    set(axspec,'Ylim',[10*log10(threshold*0.01) 0]);
  276.    set(axspec,'Xlim', [fmin*fs fmax*fs]);
  277.   end
  278.   xlabel('Energy spectral density');
  279.   Nsp=length(spec);
  280.   set(gca,['X' TickLabelStr],[],'view',[-90 90]);

  281. end;

  282.   if iscolorbar,                     % Is there a color bar ?
  283.    PositionTfr=get(axtfr,'Position');
  284.    set(axtfr,'Position',PositionTfr-[0 0 0.03 0]);
  285.    set(axcb, 'Position',[PositionTfr(1)+PositionTfr(3)-0.01,...
  286.                          PositionTfr(2) 0.01 PositionTfr(4)]);
  287.    axes(axcb);
  288.    Ncolors=length(colormap);
  289.    [cmin,cmax]=caxis; colorvect=linspace(cmax,cmin,Ncolors);
  290.    imagesc(colorvect'); axis('off');
  291.    if issig,                          % there is a signal
  292.     PositionSig=get(axsig,'Position');
  293.     set(axsig,'Position',PositionSig-[0 0 0.03 0]);
  294.    end
  295.   end

  296. axes(axtfr);                         % Display the tfr
  297. if (tcol==1),
  298.   plot(freqr,tfr(1:Nf2));
  299.   set(axtfr,'Xlim', [fmin*fs fmax*fs]);

  300. else
  301.   if strcmp(computer,'MAC'),
  302.    tfr=flipud(tfr(1:Nf2,:));
  303.   else
  304.    tfr=tfr(1:Nf2,:);
  305.   end;

  306.   indmin=find(tfr<mini);
  307.   tfr(indmin)=mini*ones(1,length(indmin));

  308.   indmax=find(tfr>maxi);
  309.   tfr(indmax)=maxi*ones(1,length(indmax));
  310.   if isaffine & (display==2),
  311.    fprintf('imagesc does not support non-linear scales for axes. Replaced by pcolor.\n');
  312.    display=3;
  313.   end;

  314.   if display==1,                  % contour
  315.    contour(ts,freqr,tfr,levels);  % contour(tfr,levels,ts,freqr);
  316.    set(axtfr,'Ylim', [fmin*fs fmax*fs]);
  317.    DisplayStr=', contour';
  318.   elseif display==2,              % imagesc
  319.    if linlogtfr==0,
  320.     imagesc(ts,freqr,tfr); axis('xy');
  321.    else
  322.     imagesc(ts,freqr,log10(tfr));axis('xy');
  323.    end
  324.    set(axtfr,'Ylim', [fmin*fs fmax*fs]);
  325.    DisplayStr=', imagesc';

  326.   elseif display==3,              % pcolor
  327.    if linlogtfr==0,
  328.     pcolor(ts,freqr,tfr); shading interp;
  329.    else
  330.     pcolor(ts,freqr,log10(tfr)); shading interp;
  331.    end
  332.    set(axtfr,'Ylim', [fmin*fs fmax*fs]);
  333.    DisplayStr=', pcolor';

  334.   elseif display==4,              % surf
  335.    if linlogtfr==0,
  336.     surf(ts,freqr,tfr); shading interp;
  337.     zlabel('Amplitude');
  338.     axis([ts(1) ts(tcol) fmin*fs fmax*fs mini maxi]);
  339.    else
  340.     surf(ts,freqr,log10(tfr)); shading interp;
  341.     zlabel('Positive values');
  342.     axis([ts(1) ts(tcol) fmin fmax log10(mini) log10(maxi)]);
  343.    end
  344.    DisplayStr=', surf';

  345.   elseif display==5,              % mesh
  346.    if linlogtfr==0,
  347.     mesh(ts,freqr,tfr); shading interp;
  348.     zlabel('Amplitude');
  349.     axis([ts(1) ts(tcol) fmin*fs fmax*fs mini maxi]);
  350.    else
  351.     mesh(ts,freqr,log10(tfr));shading interp;
  352.     zlabel('Positive values');
  353.     axis([ts(1) ts(tcol) fmin*fs fmax*fs log10(mini) log10(maxi)]);
  354.    end
  355.    DisplayStr=', mesh';
  356.   end

  357. % Define the title and check the input arguments depending on 'method'

  358. method=method(4:length(method));

  359. if nargin==5, % if there is no additional parameters, do the best.
  360.   title([method, LinLogStr,DisplayStr,...
  361.          ', Threshold=',num2str(threshold),'%']);

  362. elseif strcmp(method,'WV'  ) | strcmp(method,'MH') | ...
  363.         strcmp(method,'PAGE'), % no parameters
  364.   title([method,', Nf=',num2str(Nf2), LinLogStr,DisplayStr,...
  365.         ', Threshold=',num2str(threshold),'%']);

  366. elseif strcmp(method,'PWV'  )|strcmp(method,'PMH'  )| ...
  367.         strcmp(method,'SP'   )|strcmp(method,'PPAGE')| ...
  368.         strcmp(method,'RSP'  )|strcmp(method,'RPPAG')| ...
  369.         strcmp(method,'RPWV' )|strcmp(method,'RPMH' ),
  370.   h=p1;[hrow,hcol]=size(h); Lh=(hrow-1)/2; % one parameter
  371.   if (hcol~=1)|(rem(hrow,2)==0),
  372.    error('h must be a smoothing window with odd length'); end;
  373.   title([method, ', Lh=',num2str(Lh), ', Nf=',num2str(Nf2),...
  374.         LinLogStr, DisplayStr,', Threshold=',num2str(threshold),'%']);

  375. elseif strcmp(method,'STFT'), % short-time fourier transform case
  376.   h=p1;[hrow,hcol]=size(h); Lh=(hrow-1)/2;
  377.   if (hcol~=1)|(rem(hrow,2)==0),
  378.    error('h must be a smoothing window with odd length'); end;
  379.   title(['|',method, '|^2, Lh=',num2str(Lh),...
  380.         ', Nf=',num2str(Nf2), LinLogStr, DisplayStr,', Thld=',...
  381.         num2str(threshold),'%']);

  382. elseif strcmp(method,'SPWV' ) | strcmp(method,'MHS'  )| ...
  383.         strcmp(method,'RSPWV') | strcmp(method,'RMHS' )| ...
  384.         strcmp(method,'ZAM'  ) | strcmp(method,'RIDBN')|...
  385.         strcmp(method,'BJ'   ) | strcmp(method,'RIDB' )| ...
  386.         strcmp(method,'RIDH' ) | strcmp(method,'RIDT' ),
  387.   g=p1; [grow,gcol]=size(g); Lg=(grow-1)/2;
  388.   if (gcol~=1)|(rem(grow,2)==0),
  389.    error('g must be a smoothing window with odd length'); end;
  390.   h=p2; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
  391.   if (hcol~=1)|(rem(hrow,2)==0),
  392.    error('h must be a smoothing window with odd length'); end;
  393.   title([method,', Lg=',num2str(Lg),', Lh=',num2str(Lh),...
  394.          ', Nf=',num2str(Nf2),LinLogStr, DisplayStr,...
  395.          ', Threshold=',num2str(threshold),'%']);

  396. elseif strcmp(method,'MMCE'),
  397.   h=p1;[hrow,hcol]=size(h); Lh=(hrow-1)/2;
  398.   if (rem(hrow,2)==0),
  399.    error('h must be a smoothing window with odd length'); end;
  400.   title([method, ', Lh=',num2str(Lh), ', Nf=',num2str(Nf2),...
  401.         LinLogStr, DisplayStr,', Threshold=',num2str(threshold),'%']);

  402. elseif strcmp(method,'CW' ) | strcmp(method,'BUD'),
  403.   g=p1; [grow,gcol]=size(g); Lg=(grow-1)/2;
  404.   if (gcol~=1)|(rem(grow,2)==0),
  405.    error('g must be a smoothing window with odd length'); end;
  406.   h=p2; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
  407.   if (hcol~=1)|(rem(hrow,2)==0),
  408.    error('h must be a smoothing window with odd length'); end;
  409.   sigma=p3;
  410.   title([method,', Lg=',num2str(Lg),', Lh=',num2str(Lh),...
  411.          ' sigma=',num2str(sigma),', Nf=',num2str(Nf2),...
  412.          LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);

  413. elseif strcmp(method,'GRD')
  414.   g=p1; [grow,gcol]=size(g); Lg=(grow-1)/2;
  415.   if (gcol~=1)|(rem(grow,2)==0),
  416.    error('g must be a smoothing window with odd length'); end;
  417.   h=p2; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
  418.   if (hcol~=1)|(rem(hrow,2)==0),
  419.    error('h must be a smoothing window with odd length'); end;
  420.   title([method,', Lg=',num2str(Lg),', Lh=',num2str(Lh),...
  421.          ', rs =',num2str(p3), ', M/N =',num2str(p4), ...
  422.          ', Nf =',num2str(Nf2), ...
  423.          LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);

  424. elseif strcmp(method,'MSC' ) | strcmp(method,'RMSC' )
  425.   f0T=p1; if (f0T<=0), error('f0T must be positive'); end;
  426.   title([method, ', f0T=',num2str(f0T), ', Nf=',num2str(Nf2),...
  427.         LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);

  428. elseif strcmp(method,'RGAB' )
  429.   Nh=p1; if (Nh<=0), error('Nh must be positive'); end;
  430.   title([method, ', Nh=',num2str(Nh), ', Nf=',num2str(Nf2),...
  431.         LinLogStr, DisplayStr, ', Threshold=',num2str(threshold),'%']);


  432. elseif strcmp(method,'DFLA' ) | strcmp(method,'UNTER' )| ...
  433.         strcmp(method,'BERT' ),
  434.   N=p1;
  435.   if (N<=0),                error('N must be positive'); end;
  436.   title([method, ', N=',num2str(N), LinLogStr, DisplayStr, ', Threshold=',...
  437.          num2str(threshold), '%']);  

  438. elseif strcmp(method,'SCALO'),
  439.   Nh0=p1; N=p2;
  440.   if (Nh0<0),               error('Nh0 must be positive'); end;
  441.   if (N<=0),                error('N must be positive'); end;
  442.   if (Nh0>0),
  443.    title([method, ', Morlet wavelet, Nh0=', num2str(Nh0), ...
  444.          ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
  445.          num2str(threshold), '%']);  
  446.   else
  447.    title([method, ', Mexican hat, N=',num2str(N), LinLogStr, DisplayStr, ...
  448.           ', Thld=', num2str(threshold), '%']);  
  449.   end

  450. elseif strcmp(method,'ASPW'),
  451.   Nh0=p1; Ng0=p2; N=p3;
  452.   if (Nh0<0),               error('Nh0 must be positive'); end;
  453.   if (Ng0<0),               error('Ng0 must be positive'); end;
  454.   if (N<=0),                error('N must be positive'); end;
  455.   if (Nh0>0),
  456.    title([method, ', Morlet wlt, Nh0=', num2str(Nh0), ', Ng0=',...
  457.          num2str(Ng0), ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
  458.          num2str(threshold), '%']);  
  459.   else
  460.    title([method, ', Mexican hat, Ng0=',num2str(Ng0),...  
  461.          ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
  462.          num2str(threshold), '%']);  
  463.   end

  464. elseif strcmp(method,'SPAW'),
  465.   K=p1; Nh0=p2; Ng0=p3; N=p4;
  466.   if (Nh0<0),               error('Nh0 must be positive'); end;
  467.   if (Ng0<0),               error('Ng0 must be positive'); end;
  468.   if (N<=0),                error('N must be positive'); end;
  469.   if (Nh0>0),
  470.    title([method, ', K=', num2str(K), ', Morlet wlt, Nh0=',...
  471.          num2str(Nh0), ', Ng0=',...
  472.          num2str(Ng0), ', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
  473.          num2str(threshold), '%']);  
  474.   else
  475.    title([method, ', K=', num2str(K), ', Mexican hat, Ng0=',...
  476.          num2str(Ng0),', N=',num2str(N), LinLogStr, DisplayStr, ', Thld=',...
  477.          num2str(threshold), '%']);  
  478.   end


  479. elseif strcmp(method,'GABOR'),
  480.   N=p1; Q=p2; h=p3; [hrow,hcol]=size(h); Lh=(hrow-1)/2;
  481.   if (hcol~=1)|(rem(hrow,2)==0),
  482.    error('h must be a smoothing window with odd length'); end;
  483.   title([method, ', Lh=',num2str(Lh), ', Nf=',...
  484.          num2str(Nf2),', N=',num2str(N),', Q=',num2str(Q),...
  485.          LinLogStr, DisplayStr, ', Thld=',num2str(threshold),'%']);

  486. end;
  487. end

  488. % add the correct legend on the axes
  489. if unitHz==1,
  490. xlabel('Time [s]'); ylabel('Frequency [Hz]');
  491. elseif unitHz==2,
  492. xlabel('Time [ms]'); ylabel('Frequency [kHz]');
  493. elseif unitHz==3,
  494. xlabel('Time [祍]'); ylabel('Frequency [MHz]');
  495. end


  496. if (isgridsig & issig),        % Updating of the grids
  497. axes(axsig); grid on   
  498. elseif (~isgridsig & issig),
  499. axes(axsig); grid off
  500. end
  501. if (isgridspec & isspec),
  502. axes(axspec); grid on
  503. elseif (~isgridspec & isspec),
  504. axes(axspec); grid off
  505. end

  506. if (isgridtfr),            % upating of the grid on the tfr
  507. axes(axtfr); grid on  
  508. elseif (~isgridtfr),
  509. axes(axtfr); grid off
  510. end
复制代码

评分

1

查看全部评分

发表于 2007-4-28 20:51 | 显示全部楼层

回复 #3 zhangnan3509 的帖子

就是我在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合法的输入,还请大侠不吝赐教。
发表于 2008-1-14 15:42 | 显示全部楼层
有谁知道?感谢ing~~~~~~~~~~
发表于 2010-9-9 09:47 | 显示全部楼层
也有同样的疑问
发表于 2010-9-9 09:47 | 显示全部楼层
多学习学习?小波方差是怎么做的?

点评

请善用编辑功能!  发表于 2010-9-9 11:27
发表于 2010-9-9 14:30 | 显示全部楼层
回复 jasonwu 的帖子
建议你在版内搜索“连续小波变换”

   
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-26 01:36 , Processed in 0.058517 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表