Manual spectrogram creation without using spectrogram command
显示 更早的评论
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
采纳的回答
Mathieu NOE
2021-1-18
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 个评论
Roozbeh
2023-10-11
It seems that "Decimate" is available only in "Signal_Processing_Toolbox"!
Mathieu NOE
2023-10-11
yes , remove that portion of the code if needed
but it's worth having the Signal Processing Toolbox if you are in that job !
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Time-Frequency Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
