How to improve spectrogram plot?

Hello everyone,

I have a general question about spectogram. I have done a plot of spetrogram referent of a strong earthquake. The goal is to try to identify some other event occurre seconds later which overlapp the signal of the earthquake. I used the following command to plot:

[Stot,freqcy,temps] = spectrogram(signal_n1_spect_ev1, 300, 299, 1000, Fs_ev1);

I saw that after plot, the plot looks like to have many small interferences (i believe), in which the picture not stay very clear and solid. I have tried other values of parameters, but that was my best. Does someone have any suggestion to improve it?

Thanks in advance.

Guilherme

[Mathieu NOE]

hello

it would be more efficient to share the code + some data

[Guilherme de Melo]

Hi,

Thanks. I converted the "signal_n1_spect_ev1" variable to txt. The original data is a seismological format and you would need to install some extensions to read and extract the informations. You can find attached the signal and another txt with the timeaxis. I did like this to plot:

Fs_ev1=100;

[Stot,freqcy,temps] = spectrogram(signal_n1_spect_ev1, 300, 299, 1000, Fs_ev1);

%figure(2)
%plot spectogram
subplot(2,1,2);
%time_days = start_time + temps ./ (24*60*60);
imagesc(timeaxis1_ev1,freqcy,10*log10(abs(Stot)));
datetick('x','HH:MM:SS');
set(gca(gcf),'ydir','normal');
caxis([str2num(spec_color_scale1)]);
axis tight
xlabel('Time');
ylabel('Frequency, Hz');
colormap(parula);

You can find attached also my current version of the spectrogram and signal.

Answer 1

[Mathieu NOE]

I played a bit with some parameters and changed from imagesc to surf (pcolor would also work) to allow for optionnal frequency display in log scale (if that helps) : imagesc() Y Axis Log Scale Not Working (Help!) - MATLAB Answers - MATLAB Central

I don't know which details you are looking at , but you can always restrict the x axis to where you want and adapt also the min / max values used in caxis to have best rendering

some result obtained so far :

code :

NB : i had to change the variables names as this is how they appeared when I loaded the txt files - simply put back your original VN

Fs_ev1=100;
nfft = 512*2;
noverlap = round(0.95*nfft);
win = kaiser(nfft,18);

[Stot,freqcy,temps] = spectrogram(signaln1spectev1, win, noverlap, nfft, Fs_ev1);
%figure(2)
%plot spectogram
% subplot(2,1,2);
%time_days = start_time + temps ./ (24*60*60);
% imagesc(timeaxis,freqcy,10*log10(abs(Stot)));
% datetick('x','HH:MM:SS');
% set(gca(gcf),'ydir','normal');

s = surf(temps,freqcy,10*log10(abs(Stot)));
set(s,'edgecolor','none');
datetick('x','HH:MM:SS');
view([0 90])
set(gca, 'YScale', 'log');
colorbar('vert');
caxis([40 65])
ylim([0.5 Fs_ev1/2])

xlabel('Time');
ylabel('Frequency, Hz');
colormap(jet);
shading interp

[Mathieu NOE]

for non stationnary data , I would also recommend to try : cwt

there 's an example based on earthquake signals

[Guilherme de Melo]

Nice, it is very better to see some weak signals. Thanks for the example. Apparently, the datetick is not working with the plot using surf. :-(

[Guilherme de Melo]

Do you whay it keep that clear space at the end of the plot using the log Y-axis?

[Mathieu NOE]

workaround for datetick :

Fs_ev1=100;
nfft = 512*2;
noverlap = round(0.95*nfft);
win = kaiser(nfft,18);

[Stot,freqcy,temps] = spectrogram(signaln1spectev1, win, noverlap, nfft, Fs_ev1);
%figure(2)
%plot spectogram
% subplot(2,1,2);
%time_days = start_time + temps ./ (24*60*60);
% imagesc(timeaxis,freqcy,10*log10(abs(Stot)));
% datetick('x','HH:MM:SS');
% set(gca(gcf),'ydir','normal');

temps = temps - temps(1);
td = datevec(temps/(24*60*60)); % Convert seconds to serial date number
td = datetime(td);
td.Format = 'HH:mm:ss';

s = pcolor((1:numel(temps)),freqcy,10*log10(abs(Stot)));
set(s,'edgecolor','none');
set(gca,'TickDir','out'); % The only other option is 'in'
nticks = 4; % select how many xtick labels to display
ind = round(linspace(1,numel(temps),nticks));
XTickLabels = string(td(ind,1))';
set(gca, 'XTick', ind, 'XTickLabel', XTickLabels);

view([0 90])
set(gca, 'YScale', 'log');
colorbar('vert');
caxis([40 65])
ylim([0.5 Fs_ev1/2])

xlabel('Time');
ylabel('Frequency, Hz');
colormap(jet);
shading interp

[Mathieu NOE]

hello again

problem solved ?

do you mind accepting my answer ? tx

[Guilherme de Melo]

Hi Mathieu, yes, your helped with my issue. Thanks for accept.

[Mathieu NOE]

as always, my pleasure !