how can I convert linear figure's axis to logarithmic?

Hi,

I have following Figure#1, I want to convert only y-axis to logaritmic scale. I tried several ways but each time outcome is not normal. Can you please tell me how to achive y-axis of Figure#2. Thanks for you help.

% this my figure script. 
figure;
imagesc(t, fliplr(f), ST_normalized)

previously I used the follwing command, it did not work. I want to have the image of Figure#2 in y-axis.

for instance, "10" in Figure#1 shoud be "10^1" in a new scaled figure.

set(gca, 'YScale', 'log')

[Mathieu NOE]

imagesc() Y Axis Log Scale Not Working (Help!) - MATLAB Answers - MATLAB Central

Answer 1

[Star Strider]

I am not certain what the problem is, however you can change the ruler properties with the Exponent property.

N = 100;
x = linspace(0, 1, N);
y = linspace(1, 50, N);

figure
plot(x, y)
grid

figure
plot(x, y)
grid
Ax = gca;
Ax.YScale = 'log';
Ax.YAxis.Exponent = 1;

.

[Adnan]

Thanks for that, @Star Strider. Can you please make the same scaling for the following command? I have attached the first figure (for the linear y-axis); after scaling in the log, some of the data is missing. Blow 1 and above data are missing. How can I fix that? Thanks again.

imagesc(.)

[Star Strider]

I don’t have that file, so I can’t work with the imagesc version of it, and so can’t determine what the problem is. I might be able to help with that if I have the data and the code you use to plot it.

The best I can do otherwise is my previous example.

EDIT — (26 Mar 2025 at 14:20)

The best optiom may be to abandon imagesc and use surf instead. You can do this with the MATLAB spectrogram functions (spectrogram, pspectrum) by getting the outputs.

Using pspectrum —

LD = load('handel.mat')

LD = struct with fields:
     y: [73113x1 double]
    Fs: 8192

signal = LD.y;
Fs = LD.Fs;

[p,f,t] = pspectrum(signal, Fs, 'spectrogram');

figure
surf(t,f,p, EdgeColor='none')
colormap(turbo)
colorbar
xlabel('Time (s)')
ylabel('Frequency (Hz)')
view(0,90)
[x0,x1] = bounds(t);
[y0,y1] = bounds(f);
axis([x0 x1  y0 y1])
Ax = gca;
Ax.YScale = 'log';
Ax.YAxis.Exponent = 1;

●

This may be the better option.

Meanwhile, I cannot address the specific problems with your data without having them to work with.

.

[Adnan]

I am grateful for your support @Star Strider. If I share the code, would you be able to examine the plot section concerning my question? I have dedicated some time to this issue but have not found a solution for applying a logarithmic scale to the y-axis. I believe this might be linked to the "imagesc()" command.

[Star Strider]

My pleasure!

Yes. I will see what I can do. I will probably use a version of my code with your data.

[Adnan]

Thanks again @Star Strider, please see it attached here.

[Star Strider]

My pleasure!

Is this getting closer to what you want?

uz = unzip('matlab.zip')

uz = 1x3 cell array
    {'matlab/Ertq.txt'}    {'matlab/Main_run.m'}    {'matlab/s_trans.m'}

%% Main RUN 
%% 
%% 
% % clear all 
% % clc
% load Ertq.txt
load(uz{1})
whos

  Name         Size            Bytes  Class     Attributes

  Ertq      8372x1             66976  double              
  ans          1x34               68  char                
  uz           1x3               456  cell                

EQ=Ertq(:,1);
dt=0.005;
%% 
minfreq=0;           
maxfreq=length(EQ) ;   
samplingrate=dt ;    
[ST,t,f]=s_trans(EQ,minfreq,maxfreq,samplingrate);

Maxfreq < 0 or > Nyquist. Setting maxfreq = 4186 a/nMinfreq = 0 a/nMaxfreq = 4186 a/nSampling Rate (time   domain) = 5.000000e-03 a/nSampling Rate (freq.  domain) = 1 a/nThe length of the timeseries is 8372 points
The number of frequency voices is 4187
Calculating analytic signal (using Hilbert transform)
Estimated time is 22.735410Calculating S transform...
Finished Calculation

%% 
% Find the maximum amplitude at each time step (each column)
max_amplitude = max(abs(ST), [], 1); % Max value along each column
%% 
% Normalize by dividing each spectral amplitude by the highest amplitude at that time step
ST_normalized = abs(ST)./ max_amplitude;
%% 
% ST_normalized = abs(ST)./1;
ST_normalized=flipud(ST_normalized);
figure;
% imagesc(t,f, ST_normalized)
%% 
% % imagesc(t, fliplr(f), ST_normalized)
surf(t, fliplr(f), ST_normalized, EdgeColor='none')
ax = gca;
ax.YDir = 'normal';
ylabel('Frequency, Hz')
xlabel('Time (s)')
colorbar;
colormap(jet(256));
view(0,90)
axis('tight')
ax.YScale = 'log';
ax.YAxis.Exponent = 1;

[fmin,fmax] = bounds(f)                             % Information: ‘f’ Limits (Added, Can Be Deleted)

fmin = 0

fmax = 100

% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


function [st,t,f] = s_trans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate)
  verbose = true;           
  removeedge= false; 
  analytic_signal= isreal(timeseries); 
  factor = 3; 
  % 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(newline), 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 
  % Make sure it is a 1-dimensional array 
  assert(numel(timeseries)>1,'Please enter a *vector* of data, not a scalar')
  assert(isvector(timeseries),'Please enter a *vector* of data, not matrix') 
  % Ensure is column vector 
  if ~iscolumn(timeseries), timeseries=timeseries.';	end 
  % use defaults for input variables not provided
  switch nargin
    case 1 
      minfreq = 0; 
      maxfreq = fix(length(timeseries)/2); 
      samplingrate=1; 
      freqsamplingrate=1; 
    case 2 
      maxfreq = fix(length(timeseries)/2); 
      samplingrate=1; 
      freqsamplingrate=1; 
    case 3  
      samplingrate=1; 
      freqsamplingrate=1; 
      %[minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries); 
    case 4    
      freqsamplingrate=1; 
      %[minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries); 
    case 5 
      %[minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries); 
    otherwise
      if verbose, disp('Error in input arguments: using defaults'), end
      minfreq = 0; 
      maxfreq = fix(length(timeseries)/2); 
      samplingrate=1; 
      freqsamplingrate=1; 
  end
  [minfreq,maxfreq,samplingrate,freqsamplingrate]= check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries); 
  if verbose  
    fprintf('Minfreq = %d a/n',minfreq)
    fprintf('Maxfreq = %d a/n',maxfreq) 
    fprintf('Sampling Rate (time   domain) = %d a/n',samplingrate)
    fprintf('Sampling Rate (freq.  domain) = %d a/n',freqsamplingrate)
    fprintf('The length of the timeseries is %d points\n',length(timeseries))
  end %END OF INPUT VARIABLE CHECK 
  % 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 + (0:spe_nelements-1)*freqsamplingrate)/(samplingrate*length(timeseries)); 
  if verbose, fprintf('The number of frequency voices is %d\n',spe_nelements),end 
  % The actual S Transform function is here: 
  st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor);  
  % if no output, then plot amplitude spectrum 
  if nargout==0  
    if verbose, disp('Plotting pseudocolor image'),end 
    imagesc(t,f,abs(st))
    colormap('jet')
    colorbar
  end 
end 
 
function st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor)
n=length(timeseries); 
original = timeseries; 
if removeedge 
  if verbose, disp('Removing trend with quadratic polynomial fit'),end 
	ind = (0:n-1).'; 
  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&[1:sh_len]; 
  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 = [1; 2*ones(fix((n-1)/2),1); ones(1-rem(n,2),1); zeros(fix((n-1)/2),1)]; 
   ts_spe(:) = ts_spe.*h(:); 
   timeseries = ifft(ts_spe); 
end   
 
% Compute FFT's 
tic;vector_fft=fft(timeseries);tim_est=toc; 
vector_fft=[vector_fft,vector_fft]; 
tim_est = tim_est*ceil((maxfreq - minfreq+1)/freqsamplingrate)   ; 
if verbose, fprintf('Estimated time is %f',tim_est),end 
 st=zeros(ceil((maxfreq - minfreq+1)/freqsamplingrate),n); 
if verbose, disp('Calculating S transform...'),end 
if minfreq == 0 
   st(1,:) = mean(timeseries)*(1&[1:1:n]); 
else 
  	st(1,:)=ifft(vector_fft(minfreq+1:minfreq+n).*g_window(n,minfreq,factor)); 
end 

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 %%% end strans function 
 
%------------------------------------------------------------------------ 
function gauss=g_window(len,freq,factor) 
 
  
  vector(1,:)=(0:len-1);
  vector(2,:)=(-len:-1);
  vector=vector.^2;     
  vector=vector*(-factor*2*pi^2/freq^2); 
  
  gauss=sum(exp(vector)); 
end 

function [minfreq,maxfreq,samplingrate,freqsamplingrate]=check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries) 

  s = size(minfreq); 
  l = max(s); 
  if l > 1   
  if verbose, disp('Array of inputs accepted.'), end 
  temp=minfreq; 
  minfreq = temp(1);
  if l > 1,  maxfreq = temp(2); end
  if l > 2,  samplingrate = temp(3); end
  if l > 3,  freqsamplingrate = 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, fprintf('Maxfreq < 0 or > Nyquist. Setting maxfreq = %d a/n',maxfreq),end 
  end 
  if minfreq>maxfreq  
    temporary = minfreq; 
    minfreq = maxfreq; 
    maxfreq = 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 
end  % check_input
 
function st_help 
  disp(' ') 
	disp('st()  HELP COMMAND') 
	disp('st() returns  - 1 or an error message if it fails') 
	disp('USAGE::    [localspectra,timevector,freqvector] = 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    

.

[Adnan]

That is perfect. Thanks a lot again @Star Strider. May I ask? your script freeze my laptop under another input record (e.g., earthqauke). It makes my laptop freeze which is not normal (Lenovo, i7, 9th Gen). I even cannot enter Google while running your script. Can you please make your scrip a bit simple so that I can run multiple times otherwise it will take ages! :).

Thanks again for your time sir.

regards

[Star Strider]

My pleasure!

I am not sure what the problem could be with your laptop,.other than its needing more memory. It could be doing memory swaps on your HDD/SSD, and that would slow it down.

Other than the unzip call to import your data, the only changes I made to your original code were to add these lines:

surf(t, fliplr(f), ST_normalized, EdgeColor='none')

view(0,90)

axis('tight')

ax.YScale = 'log';

ax.YAxis.Exponent = 1;

I did not look closely at your data, however the only option I can think of would be to reduce the number of elements in the input vectors to your ‘strans’ function, possibly using the Signal Processiing Toolbox resample function. It depends on how fine the resolution needs to be.

Another approach would be to calculate ‘t’, ‘f’, and ‘ST’ separately, save them to a .mat file, exit the code that created them, and then load the .mat file to do the plotting. This would reduce the memory requirements, since the original data would not be necessary to have in memory aloing with the calculated matrices.

Othere than that, see the documentation on the memory function (especailly the See Also section) to see if you can make the memory use a bit more efficient.

.

[Adnan]

Thanks for your help sir.

[Star Strider]

As always, my pleasure!

I learn from answering your questions.

[Adnan]

Dear @Star Strider, I would like to express my gratitude for your assistance once more. I have a question to ask. After implementing your code, my computer froze. I have attempted several methods to resolve this issue, but none have been successful. Could you please assist me? I am curious if I could send you my code for you to run on your end. Is the freezing issue related to the use of a logarithmic axis? I do not encounter any problems when using a linear axis plot. However, when I switch to a logarithmic axis, my system freezes. I would greatly appreciate your help with this MATLAB plotting issue. Thank you once again.

Regards

[Star Strider]

Again, as always, my pleasure!

I doubt the logarithmic transformation would cause the freeze. It is more likely the size of the matrices, although I cannot be certain.

I will run your code (please iinclude any data files) on MATLAB Online and see if there is a problem.

[Adnan]

I appreciate your assistance again. Please find the attached file. Conduct the simulation with both logarithmic and linear axes. The results show a notable difference in simulation times for each case. Thank you.I would greatly appreciate your assistance in resolving this issue.

[Star Strider]

MATLAB Online is clkearly having problems with your code, however is throwing no errors or warnings.

A whos call just before the second plot returns:

Name Size Bytes Class Attributes

EQ 18583x1 148664 double

KUT_090 18583x2 297328 double

ST 9292x18583 2762771776 double complex

ST_normalized 9292x18583 1381385888 double

dt 1x1 8 double

f 1x9292 74336 double

max_amplitude 1x18583 148664 double

maxfreq 1x1 8 double

minfreq 1x1 8 double

samplingrate 1x1 8 double

t 1x18583 148664 double

time 18583x1 148664 double

uz 1x3 498 cell

zipf 1x30 60 char

zipfn 1x87 174 char

Ths ‘ST’ and ‘ST_normalized’ arrays together require 4,144,157,664 bytes. That could be a problem.

I am running this script:

% % % Adnan_2025_04_08.m

clear

close all

zipfn = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1831478/MATLAB_FILES.zip';

zipf = websave('MATLAB_FILES.zip',zipfn);

uz = unzip(zipf)

% % fprintf(1,['\n' repmat('=',1,72) '\n'])

% %

% % type(uz{2})

% %

% % fprintf(1,['\n' repmat('=',1,72) '\n'])

% %

% % type(uz{3})

% %

% % fprintf(1,['\n' repmat('=',1,72) '\n'])

% return

tic

%% Main RUN

%%

% clear all;

% clc

%%

% load KUT_090.txt

load(uz{1})

EQ=KUT_090(:,2);

time=KUT_090(:,1);

dt = mean(diff(time));

%%

minfreq=0;

maxfreq=length(EQ) ;

samplingrate=dt ;

[ST,t,f]=s_trans_st(EQ,minfreq,maxfreq,samplingrate);

max_amplitude = max(abs(ST), [], 1);

ST_normalized = abs(ST)./ max_amplitude;

ST_normalized=flipud(ST_normalized);

%%

figure;

imagesc(t, fliplr(f), ST_normalized)

ax = gca;

ax.YDir = 'normal';

ylabel('Frequency, Hz')

xlabel('Time (s)')

colorbar;

colormap(jet(256));

%%

%%

%%

figure;

surf(t, fliplr(f), ST_normalized, EdgeColor='none')

ax = gca;

ax.YDir = 'normal';

ylabel('Frequency, Hz')

xlabel('Time (s)')

colorbar;

colormap(jet(256));

view(0,90)

axis('tight')

ax.YScale = 'log';

ax.YAxis.Exponent = 1;

yline(0.305, 'c--', 'LineWidth', 2); % sys_ini, FF1

yline(0.286, 'black--', 'LineWidth', 2); % sys_inc_m, FF1

xline(70.03, 'r--', 'LineWidth', 2); %sys_ini, FF1

xline(172.94, 'black--', 'LineWidth', 2) ; %sys_inc_m ,,, FF1

toc

% ========================================================================

function [st,t,f] = s_trans_st(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate)

verbose = true;

removeedge= false;

analytic_signal= isreal(timeseries);

factor = 3;

% 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(newline), 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

% Make sure it is a 1-dimensional array

assert(numel(timeseries)>1,'Please enter a *vector* of data, not a scalar')

assert(isvector(timeseries),'Please enter a *vector* of data, not matrix')

% Ensure is column vector

if ~iscolumn(timeseries), timeseries=timeseries.'; end

% use defaults for input variables not provided

switch nargin

case 1

minfreq = 0;

maxfreq = fix(length(timeseries)/2);

samplingrate=1;

freqsamplingrate=1;

case 2

maxfreq = fix(length(timeseries)/2);

samplingrate=1;

freqsamplingrate=1;

case 3

samplingrate=1;

freqsamplingrate=1;

%[minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

case 4

freqsamplingrate=1;

%[minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

case 5

%[minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

otherwise

if verbose, disp('Error in input arguments: using defaults'), end

minfreq = 0;

maxfreq = fix(length(timeseries)/2);

samplingrate=1;

freqsamplingrate=1;

end

[minfreq,maxfreq,samplingrate,freqsamplingrate]= check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

if verbose

fprintf('Minfreq = %d a/n',minfreq)

fprintf('Maxfreq = %d a/n',maxfreq)

fprintf('Sampling Rate (time domain) = %d a/n',samplingrate)

fprintf('Sampling Rate (freq. domain) = %d a/n',freqsamplingrate)

fprintf('The length of the timeseries is %d points\n',length(timeseries))

end %END OF INPUT VARIABLE CHECK

% 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 + (0:spe_nelements-1)*freqsamplingrate)/(samplingrate*length(timeseries));

if verbose, fprintf('The number of frequency voices is %d\n',spe_nelements),end

% The actual S Transform function is here:

st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor);

% if no output, then plot amplitude spectrum

if nargout==0

if verbose, disp('Plotting pseudocolor image'),end

imagesc(t,f,abs(st))

colormap('jet')

colorbar

end

end

function st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor)

n=length(timeseries);

original = timeseries;

if removeedge

if verbose, disp('Removing trend with quadratic polynomial fit'),end

ind = (0:n-1).';

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&[1:sh_len];

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 = [1; 2*ones(fix((n-1)/2),1); ones(1-rem(n,2),1); zeros(fix((n-1)/2),1)];

ts_spe(:) = ts_spe.*h(:);

timeseries = ifft(ts_spe);

end

% Compute FFT's

tic;vector_fft=fft(timeseries);tim_est=toc;

vector_fft=[vector_fft,vector_fft];

tim_est = tim_est*ceil((maxfreq - minfreq+1)/freqsamplingrate) ;

if verbose, fprintf('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&[1:1:n]);

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 %%% end strans function

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

function gauss=g_window(len,freq,factor)

vector(1,:)=(0:len-1);

vector(2,:)=(-len:-1);

vector=vector.^2;

vector=vector*(-factor*2*pi^2/freq^2);

% Compute the Gaussion window

gauss=sum(exp(vector));

end

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 > 1, maxfreq = temp(2); end

if l > 2, samplingrate = temp(3); end

if l > 3, freqsamplingrate = 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, fprintf('Maxfreq < 0 or > Nyquist. Setting maxfreq = %d a/n',maxfreq),end

end

if minfreq>maxfreq

temporary = minfreq;

minfreq = maxfreq;

maxfreq = 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

end % check_input

function st_help

disp(' ')

disp('st() HELP COMMAND')

disp('st() returns - 1 or an error message if it fails')

disp('USAGE:: [localspectra,timevector,freqvector] = 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

When I suppress the first plot by commenting it out, nothing plots and the code crashes. Setting the YScale value does not change that behaviour. The elapsed time (given by the tic and toc calls I added) is generally 20-22 seconds. Available memory could be a problem, however I would expect an error or warning in that instance. I am not certain what you are doing, so I cannot suggest any changes.

.

[Adnan]

Thanks for that. Interestingly, the whole code runs fine once I use the following figuring.

figure;
imagesc(t, fliplr(f), ST_normalized)
ax = gca;
ax.YDir = 'normal';
ylabel('Frequency, Hz')
xlabel('Time (s)')
colorbar;
colormap(jet(256));

Nonetheless, when I execute the following script, it freezes. I would like to know what is wrong with the following scripts. I think the only problem is with the following scripts. The follwing script appears to be correct, but something is causing the plotting to freeze. Is that possible we can solve this issue? Thanks again @Star Strider sir for your help.

figure;
surf(t, fliplr(f), ST_normalized, EdgeColor='none')
ax = gca;
ax.YDir = 'normal';
ylabel('Frequency, Hz')
xlabel('Time (s)')
colorbar;
colormap(jet(256));
view(0,90)
axis('tight')
ax.YScale = 'log';
ax.YAxis.Exponent = 1;