Hi I have Csv file for voltage and time data and I would like to compute spectrogram to to compute harmonics at different frequencies but my spectrogram looks so much noisy or

%%
folder = 'C:\Users\Min\Desktop\New folder (3)';
filename = '27o2.csv';
data = readtable(fullfile(folder, filename));
t = table2array(data(3:end, 1));
x = table2array(data(3:end, 2));
fs = 1 / (t(2) - t(1));
% Plot the spectrogram
figure;
spectrogram(x, 500,100,500, fs, 'yaxis');
% Customize the plot
xlabel('Time (s)');
ylabel('Frequency (MHz)');
colormap("hot")
title('Spectrogram');

colorbar;

Answer 1

[Mathieu NOE]

hello

you may play with the colormap and also with the displayed dB range to focus on your signal tones and not the background

here I prefer the jet colormap with a dB range limited to 30 dB (from the max amplitude to 30 dB lower then everything below that lower limit is displayed with dark blue color)

my code uses my own spectrogram function but you can achieve the same result with the regular spectrogram function

folder = pwd;
filename = '27o2.csv';
data = readmatrix(fullfile(folder, filename));
t = data(3:end, 1);
x = data(3:end, 2);
fs = 1/mean(diff(t));

nfft = 400; 
overlap = 0.75; 
dB_range = 30; % display this dB range on y axis

[S,F,T] = myspecgram(x, 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))) - dB_range;

% plot
figure(1); 
imagesc(T,F,sg_dB) 
caxis([min_disp_dB min_disp_dB+dB_range]);
hcb = colorbar('vert');
set(get(hcb,'Title'),'String','dB')
set(gca,'YDir','Normal') 
xlabel('Time (secs)') 
ylabel('Freq (Hz)') 
title('Specgram')
colormap('jet');
% colormap('hot');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Honey]

Thank you so much for help.

[Mathieu NOE]

my pleasure !

Answer 2

[Angelo Yeo]

Hi @Cynthia,

What do you mean by "noisy"? Your calculation gives you the correct information, i.e., sum of harmonic sinusoidals. First of all, let's take a look at how your signal looks in time domain.

%%
folder = pwd;
filename = '27o2.csv';
data = readtable(fullfile(folder, filename));
t = table2array(data(3:end, 1));
x = table2array(data(3:end, 2));
fs = 1 / (t(2) - t(1));

figure;
plot(t, x);
axis tight

It looks like the signal is quite stationary, i.e., no big change of mean and variance. And I don't see a sudden change in the signal. Hence, let's take a look at the FFT result.

%% FFT
T = 1/fs;             % Sampling period       
L = length(x);             % Length of signal
t = (0:L-1)*T;        % Time vector
Y = fft(x);

figure;
plot(fs/L*(-L/2:L/2-1),abs(fftshift(Y)),"LineWidth",3)
title("fft Spectrum in the Positive and Negative Frequencies")
xlabel("f (Hz)")
ylabel("|fft(X)|")

Just like you said, it has harmonic components.

Now, let's give a slight change to your STFT calculation, i.e., give bigger overlap so that it can have a "smoother" result.

%% Plot the spectrogram
figure;
spectrogram(x, 500, 480, [], fs, 'yaxis');
% Customize the plot
xlabel('Time (s)');
ylabel('Frequency (MHz)');
colormap("hot")
title('Spectrogram');
colorbar;

[Honey]

yes I see that.Thank you so much.