c2019 发表于 2010-9-15 16:56

求广义S变换的程序

求广义S变换的MATLAB程序,先谢啦

杨德昌 发表于 2010-9-15 22:29

{:{27}:}借你的帖子同求!

ChaChing 发表于 2010-9-15 23:56

不懂! 但google一下不是就有了:@)
http://www.cora.nwra.com/~stockwel/index.php?module=pagemaster&PAGE_user_op=view_page&PAGE_id=9
function = st(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate)
% Returns the Stockwell Transform of the timeseries.
% Code by Robert Glenn Stockwell.
% DO NOT DISTRIBUTE
% BETA TEST ONLY
% Reference is "Localization of the Complex Spectrum: The S Transform"
% from IEEE Transactions on Signal Processing, vol. 44., number 4, April 1996, pages 998-1001.
%
%-------Inputs Needed------------------------------------------------
%
%   *****All frequencies in (cycles/(time unit))!******
%        "timeseries" - vector of data to be transformed
%-------Optional Inputs ------------------------------------------------
%
%"minfreq" is the minimum frequency in the ST result(Default=0)
%"maxfreq" is the maximum frequency in the ST result (Default=Nyquist)
%"samplingrate" is the time interval between samples (Default=1)
%"freqsamplingrate" is the frequency-sampling interval you desire in the ST result (Default=1)
%Passing a negative number will give the default ex. = st(data,-1,-1,2,2)
%-------Outputs Returned------------------------------------------------
%
% st   -a complex matrix containing the Stockwell transform.
%                       The rows of STOutput are the frequencies and the
%         columns are the time values ie each column is
%         the "local spectrum" for that point in time
%t      - a vector containing the sampled times
%f      - a vector containing the sampled frequencies
%--------Additional details-----------------------
%   %There are several parameters immediately below that
%the user may change. They are:
%    if true prints out informational messages throughout the function.
% if true, removes a least squares fit parabola
%                and puts a 5% hanning taper on the edges of the time series.
%                This is usually a good idea.
%if the timeseries is real-valued
%                      this takes the analytic signal and STs it.
%                      This is almost always a good idea.
%   the width factor of the localizing gaussian
%                ie, a sinusoid of period 10 seconds has a
%                gaussian window of width factor*10 seconds.
%                I usually use factor=1, but sometimes factor = 3
%                to get better frequency resolution.
%   Copyright (c) by Bob Stockwell
%   $Revision: 1.2 $$Date: 1997/07/08$


% This is the S transform wrapper that holds default values for the function.
TRUE = 1;
FALSE = 0;
%%% DEFAULT PARAMETERS
verbose = TRUE;         
removeedge= FALSE;
analytic_signal =FALSE;
factor = 1;
%%% END of DEFAULT PARAMETERS


%%%START OF INPUT VARIABLE CHECK
% First:make sure it is a valid time_series
%         If not, return the help message

if verbose disp(' '),end% i like a line left blank

if nargin == 0
   if verbose disp('No parameters inputted.'),end
   st_help
   t=0;,st=-1;,f=0;
   return
end

% Change to column vector
if size(timeseries,2) > size(timeseries,1)
        timeseries=timeseries';       
end

% Make sure it is a 1-dimensional array
if size(timeseries,2) > 1
   error('Please enter a *vector* of data, not matrix')
        return
elseif (size(timeseries)==) == 1
        error('Please enter a *vector* of data, not a scalar')
        return
end

% use defaults for input variables

if nargin == 1
   minfreq = 0;
   maxfreq = fix(length(timeseries)/2);
   samplingrate=1;
   freqsamplingrate=1;
elseif nargin==2
   maxfreq = fix(length(timeseries)/2);
   samplingrate=1;
   freqsamplingrate=1;
   [ minfreq,maxfreq,samplingrate,freqsamplingrate] =check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
elseif nargin==3
   samplingrate=1;
   freqsamplingrate=1;
   [ minfreq,maxfreq,samplingrate,freqsamplingrate] =check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
elseif nargin==4   
   freqsamplingrate=1;
   [ minfreq,maxfreq,samplingrate,freqsamplingrate] =check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
elseif nargin == 5
      [ minfreq,maxfreq,samplingrate,freqsamplingrate] =check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
else      
   if verbose disp('Error in input arguments: using defaults'),end
   minfreq = 0;
   maxfreq = fix(length(timeseries)/2);
   samplingrate=1;
   freqsamplingrate=1;
end
if verbose
   disp(sprintf('Minfreq = %d',minfreq))
   disp(sprintf('Maxfreq = %d',maxfreq))
   disp(sprintf('Sampling Rate (time   domain) = %d',samplingrate))
   disp(sprintf('Sampling Rate (freq.domain) = %d',freqsamplingrate))
   disp(sprintf('The length of the timeseries is %d points',length(timeseries)))

   disp(' ')
end
%END OF INPUT VARIABLE CHECK

% If you want to "hardwire" minfreq & maxfreq & samplingrate & freqsamplingrate do it here

% calculate the sampled time and frequency values from the two sampling rates
t = (0:length(timeseries)-1)*samplingrate;
spe_nelements =ceil((maxfreq - minfreq+1)/freqsamplingrate)   ;
f = (minfreq + *freqsamplingrate)/(samplingrate*length(timeseries));
if verbose disp(sprintf('The number of frequency voices is %d',spe_nelements)),end


% The actual S Transform function is here:
st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor);
% this function is below, thus nicely encapsulated

%WRITE switch statement on nargout
% if 0 then plot amplitude spectrum
if nargout==0
   if verbose disp('Plotting pseudocolor image'),end
   pcolor(t,f,abs(st))
end


return


%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


function st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor);
% Returns the Stockwell Transform, STOutput, of the time-series
% Code by R.G. Stockwell.
% Reference is "Localization of the Complex Spectrum: The S Transform"
% from IEEE Transactions on Signal Processing, vol. 44., number 4,
% April 1996, pages 998-1001.
%
%-------Inputs Returned------------------------------------------------
%         - are all taken care of in the wrapper function above
%
%-------Outputs Returned------------------------------------------------
%
%        ST    -a complex matrix containing the Stockwell transform.
%                       The rows of STOutput are the frequencies and the
%                       columns are the time values
%
%
%-----------------------------------------------------------------------

% Compute the length of the data.
n=length(timeseries);
original = timeseries;
if removeedge
    if verbose disp('Removing trend with polynomial fit'),end
       ind = ';
    r = polyfit(ind,timeseries,2);
    fit = polyval(r,ind) ;
       timeseries = timeseries - fit;
    if verbose disp('Removing edges with 5% hanning taper'),end
    sh_len = floor(length(timeseries)/10);
    wn = hanning(sh_len);
    if(sh_len==0)
       sh_len=length(timeseries);
       wn = 1&;
    end
    % make sure wn is a column vector, because timeseries is
   if size(wn,2) > size(wn,1)
      wn=wn';       
   end
   
   timeseries(1:floor(sh_len/2),1) = timeseries(1:floor(sh_len/2),1).*wn(1:floor(sh_len/2),1);
        timeseries(length(timeseries)-floor(sh_len/2):n,1) = timeseries(length(timeseries)-floor(sh_len/2):n,1).*wn(sh_len-floor(sh_len/2):sh_len,1);

end

% If vector is real, do the analytic signal

if analytic_signal
   if verbose disp('Calculating analytic signal (using Hilbert transform)'),end
   % this version of the hilbert transform is different than hilbert.m
   %This is correct!
   ts_spe = fft(real(timeseries));
   h = ;
   ts_spe(:) = ts_spe.*h(:);
   timeseries = ifft(ts_spe);
end

% Compute FFT's
tic;vector_fft=fft(timeseries);tim_est=toc;
vector_fft=;
tim_est = tim_est*ceil((maxfreq - minfreq+1)/freqsamplingrate)   ;
if verbose disp(sprintf('Estimated time is %f',tim_est)),end

% Preallocate the STOutput matrix
st=zeros(ceil((maxfreq - minfreq+1)/freqsamplingrate),n);
% Compute the mean
% Compute S-transform value for 1 ... ceil(n/2+1)-1 frequency points
if verbose disp('Calculating S transform...'),end
if minfreq == 0
   st(1,:) = mean(timeseries)*(1&);
else
        st(1,:)=ifft(vector_fft(minfreq+1:minfreq+n).*g_window(n,minfreq,factor));
end

%the actual calculation of the ST
% Start loop to increment the frequency point
for banana=freqsamplingrate:freqsamplingrate:(maxfreq-minfreq)
   st(banana/freqsamplingrate+1,:)=ifft(vector_fft(minfreq+banana+1:minfreq+banana+n).*g_window(n,minfreq+banana,factor));
end   % a fruit loop!   aaaaa ha ha ha ha ha ha ha ha ha ha
% End loop to increment the frequency point
if verbose disp('Finished Calculation'),end

%%% end strans function

%------------------------------------------------------------------------
function gauss=g_window(length,freq,factor)

% Function to compute the Gaussion window for
% function Stransform. g_window is used by function
% Stransform. Programmed by Eric Tittley
%
%-----Inputs Needed--------------------------
%
%        length-the length of the Gaussian window
%
%        freq-the frequency at which to evaluate
%                  the window.
%        factor- the window-width factor
%
%-----Outputs Returned--------------------------
%
%        gauss-The Gaussian window
%

vector(1,:)=;
vector(2,:)=[-length:-1];
vector=vector.^2;   
vector=vector*(-factor*2*pi^2/freq^2);
% Compute the Gaussion window
gauss=sum(exp(vector));

%-----------------------------------------------------------------------

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%
function [ minfreq,maxfreq,samplingrate,freqsamplingrate] =check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries)
% this checks numbers, and replaces them with defaults if invalid

% if the parameters are passed as an array, put them into the appropriate variables
s = size(minfreq);
l = max(s);
if l > 1
   if verbose disp('Array of inputs accepted.'),end
   temp=minfreq;
   minfreq = temp(1);;
   if l > 1maxfreq = temp(2);,end;
   if l > 2samplingrate = temp(3);,end;
   if l > 3freqsamplingrate = temp(4);,end;
   if l > 4
      if verbose disp('Ignoring extra input parameters.'),end
   end;

end      
   
   if minfreq < 0 | minfreq > fix(length(timeseries)/2);
      minfreq = 0;
      if verbose disp('Minfreq < 0 or > Nyquist. Setting minfreq = 0.'),end
   end
   if maxfreq > length(timeseries)/2| maxfreq < 0
      maxfreq = fix(length(timeseries)/2);
      if verbose disp(sprintf('Maxfreq < 0 or > Nyquist. Setting maxfreq = %d',maxfreq)),end
   end
      if minfreq > maxfreq
      temporary = minfreq;
      minfreq = maxfreq;
      maxfreq = temporary;
      clear temporary;
      if verbose disp('Swapping maxfreq <=> minfreq.'),end
   end
   if samplingrate <0
      samplingrate = abs(samplingrate);
      if verbose disp('Samplingrate <0. Setting samplingrate to its absolute value.'),end
   end
   if freqsamplingrate < 0   % check 'what if freqsamplingrate > maxfreq - minfreq' case
      freqsamplingrate = abs(freqsamplingrate);
      if verbose disp('Frequency Samplingrate negative, taking absolute value'),end
   end

% bloody odd how you don't end a function

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%
function st_help
   disp(' ')
        disp('st()HELP COMMAND')
        disp('st() returns- 1 or an error message if it fails')
        disp('USAGE::    = st(timeseries)')
        disp('NOTE::   The function st() sets default parameters then calls the function strans()')
   disp(' ')
   disp('You can call strans() directly and pass the following parameters')
   disp(' **** Warning!These inputs are not checked if strans() is called directly!! ****')
        disp('USAGE::localspectra = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor) ')
   
   disp(' ')
   disp('Default parameters (available in st.m)')
        disp('VERBOSE          - prints out informational messages throughout the function.')
        disp('REMOVEEDGE       - removes the edge with a 5% taper, and takes')
   disp('FACTOR         -the width factor of the localizing gaussian')
   disp('                  ie, a sinusoid of period 10 seconds has a ')
   disp('                  gaussian window of width factor*10 seconds.')
   disp('                  I usually use factor=1, but sometimes factor = 3')
   disp('                  to get better frequency resolution.')
   disp(' ')
   disp('Default input variables')
   disp('MINFREQ         - the lowest frequency in the ST result(Default=0)')
   disp('MAXFREQ         - the highest frequency in the ST result (Default=nyquist')
   disp('SAMPLINGRATE      - the time interval between successive data points (Default = 1)')
   disp('FREQSAMPLINGRATE- the number of frequencies between samples in the ST results')
       
% end of st_help procedure   

% the s-tranform test script
len = 128;
freq = 5;   
t = 0:len-1;

% CREATE CROSS CHIRP TIME SERIES
cross_chirp = cos(2*pi*(10+t/7).*t/len) + cos(2*pi*(len/2.8-t/6.0).*t/len);
% CREATE MODULATED SIN FUNCTION TIME SERIES
mod_freq=4*cos(2*pi*t/len)+len/5;
sin_of_sin = cos(2*pi*mod_freq.*t/len);
% CREATE the Distinct sinusoids example
midfreq = 20;
lowfreq = 5;
highfreq = 45;
distinct = cos(2*pi*midfreq.*t/len);
distinct(1:len/2) = cos(2*pi*lowfreq.*t(1:len/2)/len);
distinct(20:30) = cos(2*pi*highfreq .*t(20:30)/len);
% CREATE the chirp example
chirp = cos(2*pi*(10+t/7).*t/len);



= st(sin_of_sin);
= st(chirp);
= st(cross_chirp);
= st(distinct);

contourf(st_times,st_frequencies,abs(st_matrix_chirps));
%mesh(abs(st_matrix_chirp));

c2019 发表于 2010-9-16 09:50

谢谢,
请问运行时提示function = st(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate)
   
Error: Function definitions are not permitted at the prompt or in scripts.
是什么意思,已安装时频分析软件包,还要安装其他软件包吗?


   

Happy99 发表于 2010-9-16 12:38

回复 c2019 的帖子
LS直接复制至命令窗? 学下函数使用:@)

c2019 发表于 2010-9-18 10:08

谢谢
   

zxqzxq3077 发表于 2012-4-10 20:51

求解广义S变换扣扣1733726172

石头王石头 发表于 2012-10-31 11:39

c2019 发表于 2010-9-16 09:50 static/image/common/back.gif
谢谢,
请问运行时提示function = st(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrat ...

上面的变换是广义S变换??、还是S变换?我只想用S变换输出频率和时间的曲线,谢谢
页: [1]
查看完整版本: 求广义S变换的程序