How can I plot frequency vs magnitude of the scalogram?. The input file is the M*N matrix that is a scalogram. Any one help me by providing MATLAB code.

[sharmin sathi]

@Walter Roberson sir, could you help me?.

Answer 1

[Mathieu NOE]

hello

the spectrum is either the mean or the peak value along the x direction of S

see code example below

you can change my home made spectrogram function to the already available functions (specgram or spectrogram) if you prefer ; that was for someone else wanting to see how a spectrogram is made

see the RHS plot made as the time averaged of the dB spectrogram

fs=2000; % Sampling rate
N=5000; % Number of data points
t=[0:1:N-1]/fs; 
y=sin(250*2*pi*t)+chirp(t,100,2,500,'quadratic');

nfft = 200; 
overlap = 0.75; 
spectrogram_dB_scale = 80; % display this dB range on y axis

[S,F,T] = myspecgram(y, fs, nfft, overlap); % overlap is 75% of nfft here
sg_dB = 20*log10(abs(S)); % expressed now in dB

% saturation of the dB range :  the lowest displayed level is spectrogram_dB_scale dB below the max level
min_disp_dB = round(max(max(sg_dB))) - spectrogram_dB_scale;
sg_dB(sg_dB<min_disp_dB) = min_disp_dB;

% plot
figure(1); 
subplot(121),
imagesc(T,F,sg_dB) 
hcb = colorbar('vert');
set(get(hcb,'Title'),'String','dB')
set(gca,'YDir','Normal') 
xlabel('Time (secs)') 
ylabel('Freq (Hz)') 
title('Specgram')
colormap('jet');

subplot(122),
sdB = mean(sg_dB,2); % spectrum time averaged = mean along dim 2
plot(sdB,F) 
xlabel('Mean Amplitude (dB)') 
ylabel('Freq (Hz)') 
title('Spectrum')



function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

% main loop
fft_specgram = [];
for ci=1:spectnum
    start = (ci-1)*offset;
    sw = signal((1+start):(start+nfft)).*window;
    fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
end

% one sidded fft spectrum  % Select first half 
if rem(nfft,2)    % nfft odd
    select = (1:(nfft+1)/2)';
else
    select = (1:nfft/2+1)';
end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(nfft/2))/Fs;

end

[sharmin sathi]

Hello Mathieu NOE,

thank you so much.

I need the matlab code detail of the wavelet transform.

could you help me.

[Mathieu NOE]

hello again

I don't have the Wavelet Toolbox but I used this Fex submission

Chirplet Transform - File Exchange - MATLAB Central (mathworks.com)

you can also use this example I modified from the FEX submission :

% Test the "Chirplet_Transform" function, and compare the results of CT
% with them of STFT/CWT (just run it to see all the results mentioned above)
clear;
clc;
close all;
%% generate the simulating signal
Fs=2000; % Sampling rate
N=5000; % Number of data points
dt = 1/Fs;
t=(0:N-1)*dt; 
Sig=sin(50*2*pi*t)+chirp(t,100,2,500,'quadratic');


%% STFT
fLevel = 1024;
WinLen = 512;

[Spec,Freq] = STFT(Sig,fLevel,WinLen,Fs);

% reduction of the dB scale : only display max to min = max - Range_dB
Range_dB = 40; %  dB scale
out_dB = 20*log10(abs(Spec));
ind = out_dB< max(out_dB,[],'all')-Range_dB;
out_dB(ind) = NaN;


% plot STFT
figure
imagesc(t,Freq,out_dB);

JET = colormap(jet);
mycolor = cat(1, [1 1 1], JET(2:end,:));
colormap (JET);   % Set the current colormap Style 

box on;
colorbar off;

set(gca,'linewidth',1.5,'fontsize',15,'fontname','Times New Roman');
xlabel('Time / Sec','FontSize',15);
ylabel('Freq. / Hz','FontSize',15);
view(0,90);
set(gca,'YDir','normal');
set(gcf,'Color','w');
title('STFT');

%% CT
fLevel = 1024;
WinLen = 512;
alpha = 5; % in Hz/s
[Spec,Freq,~] = Chirplet_Transform(Sig,fLevel,WinLen,Fs,alpha);

% reduction of the dB scale : only display max to min = max - Range_dB
Range_dB = 40; %  dB scale
out_dB = 20*log10(abs(Spec));
ind = out_dB< max(out_dB,[],'all')-Range_dB;
out_dB(ind) = NaN;


% plot CT
figure
imagesc(t,Freq,out_dB);

JET = colormap(jet);
mycolor = cat(1, [1 1 1], JET(2:end,:));
colormap (JET);   % Set the current colormap Style 

box on;
colorbar off;

set(gca,'linewidth',1.5,'fontsize',15,'fontname','Times New Roman');
ylim([Freq(1),Freq(end)]);
xlabel('Time / Sec','FontSize',15);
ylabel('Freq. / Hz','FontSize',15);
view(0,90);
set(gca,'YDir','normal');
set(gcf,'Color','w');
title('CT');

%% CWT
WinLen = 128;
fLevel = 512;
[CWT, ~, Freq, Scales] = CWT_Morlet(Sig,Fs,WinLen,fLevel);

% reduction of the dB scale : only display max to min = max - Range_dB
Range_dB = 40; %  dB scale
out_dB = 20*log10(abs(CWT));
ind = (out_dB< max(out_dB,[],'all')-Range_dB);
out_dB(ind) = NaN;


% plot CWT
figure
mesh(t,Freq,out_dB);

JET = colormap(jet);
mycolor = cat(1, [1 1 1], JET(2:end,:));
colormap (JET);   % Set the current colormap Style 

box on;
colorbar off;

set(gca,'linewidth',1.5,'fontsize',15,'fontname','Times New Roman');
ylim([Freq(1),Freq(end)]);
xlabel('Time / Sec','FontSize',15);
ylabel('Freq. / Hz','FontSize',15);
view(0,90);
set(gca,'YDir','normal');
set(gcf,'Color','w');
title('CWT');

[sharmin sathi]

Hello Mathieu NOE,

thank you so much.

[Mathieu NOE]

As always, my pleasure !

[Mathieu NOE]

hello again

If my submission fullfills your request, do you mind accepting it ?

tx