Manual spectrogram creation without using spectrogram command
72 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Mohamed Medhat
le 18 Jan 2021
Commenté : Mathieu NOE
le 11 Oct 2023
I have a signal y which I need to create spectrograms of it with these parameters without using the spectrogram command
Time-domain window size: use the following values {16; 64; 256; 1024; 4096} ( meaning to split the signal 5 times, first time into windows of length 16, second into windows of length 64 and so on)
Overlap between the sliding windows: 50% (meaning if the first part is 1:16 the second is 9:24 & the third is 17:32 and so on)
FFT size: 2^13 (meaning to use 2^13 samples of the signal only)
Create a spectrogram for each window size ({16; 64; 256; 1024; 4096}) using a rectangular time-domain window.
What I'm actually stuck at is how to implement the window splitting thing, I thought of using a for loop and after splitting the windows and taking the FFTs creating the spectrogram with imagesc but not really seeing how the splitting and storing the FFT values of each window should go
0 commentaires
Réponse acceptée
Mathieu NOE
le 18 Jan 2021
hello
see below (myspecgram)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
[data,Fs] = audioread('myWAVaudiofile.wav');
channel = 1;
signal = data(:,channel);
samples = length(signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 4096; %
OVERLAP = 0.75;
% spectrogram dB scale
spectrogram_dB_scale = 100; % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 2;
if decim>1
signal = decimate(signal,decim);
Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));
[sg,fsg,tsg] = myspecgram(signal,NFFT,Fs,OVERLAP);
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% saturation of the dB range :
% saturation_dB = 60; % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([' Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function [sg,freq_vector,time] = myspecgram(signal,nfft,Fs,Overlap)
% [sg,freq_vector,time] = myspecgram(signal,NFFT,Fs,hanning(NFFT),floor(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)
dt = 1/Fs;
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
sg = [];
for i=1:spectnum
start = (i-1)*offset;
sw = signal((1+start):(start+nfft)).*window;
sg = [sg (abs(fft(sw))*4/nfft)]; % X=fft(x.*hanning(N))*4/N; % hanning only
time(i) = (start+nfft/2)*dt; % time defined as in the middle of the data buffer
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
sg = sg(select,:);
freq_vector = (select - 1)*Fs/nfft;
end
2 commentaires
Mathieu NOE
le 11 Oct 2023
yes , remove that portion of the code if needed
but it's worth having the Signal Processing Toolbox if you are in that job !
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Time-Frequency Analysis dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!