Manual spectrogram creation without using spectrogram command
43 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Mohamed Medhat
am 18 Jan. 2021
Kommentiert: Mathieu NOE
am 11 Okt. 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 Kommentare
Akzeptierte Antwort
Mathieu NOE
am 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 Kommentare
Mathieu NOE
am 11 Okt. 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 !
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Time-Frequency Analysis finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!