Search code examples
matlabwavelettime-frequency

compute S transform and its square value in matlab


let us consider following code taken from

http://www.mathworks.com/matlabcentral/fileexchange/45848-stockwell-transform--s-transform-/content/stran.m

to compute S transform, here is my code

function ST=stran(h)
% Compute S-Transform without for loops
%%% Coded by Kalyan S. Dash %%%
%%% IIT Bhubaneswar, India %%%
[~,N]=size(h); % h is a 1xN one-dimensional series
nhaf=fix(N/2);
odvn=1;
if nhaf*2==N;
    odvn=0;
end
f=[0:nhaf -nhaf+1-odvn:-1]/N;
Hft=fft(h);
%Compute all frequency domain Gaussians as one matrix
invfk=[1./f(2:nhaf+1)]';
W=2*pi*repmat(f,nhaf,1).*repmat(invfk,1,N);
G=exp((-W.^2)/2); %Gaussian in freq domain
% End of frequency domain Gaussian computation
% Compute Toeplitz matrix with the shifted fft(h)
HW=toeplitz(Hft(1:nhaf+1)',Hft);
% Exclude the first row, corresponding to zero frequency
HW=[HW(2:nhaf+1,:)];
% Compute Stockwell Transform
ST=ifft(HW.*G,[],2); %Compute voice
%Add the zero freq row
st0=mean(h)*ones(1,N);

ST=[st0;ST];

end

and consider following chirp signal

>> t = 0:0.001:2;
x = chirp(t,100,1,200,'quadratic');

i need confirmation that i am doing correctly following things

>> ST=stran(x);
>> plot(abs(ST))

?

picture is here enter image description here


Solution

  • Posting my comment as answer:

    I don't have much idea of the s-transform, but AFAIK the result of it is a 3D signal (as you can clearly see in the size of ST), so you may want to do imagesc(abs(ST)) or surf(abs(ST),'linestyle','none') instead of plot.

    In your figure you have plotted 1000 lines, that's why its so chaotic.

    Using

    imagesc(abs(ST)) enter image description here