update Matlab code

Software/Matlab STFT Program/STFT_TEST.m

@@ -0,0 +1,47 @@
+clear, clc, close all
+
+% load a .wav file
+[x, fs] = audioread('ppl2.wav');   % load an audio file
+x = x(:, 1);                        % get the first channel
+
+% define analysis parameters
+xlen = length(x);                   % length of the signal
+wlen = 512;                        % window length (recomended to be power of 2)
+hop = wlen/4;                       % hop size (recomended to be power of 2)
+nfft = 1024;                        % number of fft points (recomended to be power of 2)
+
+% perform STFT
+[S, f, t] = stft(x, wlen, hop, nfft, fs);
+
+% define the coherent amplification of the window
+K = sum(hamming(wlen, 'periodic'))/wlen;
+
+% take the amplitude of fft(x) and scale it, so not to be a
+% function of the length of the window and its coherent amplification
+S = abs(S)/wlen/K;
+
+% correction of the DC & Nyquist component
+if rem(nfft, 2)                     % odd nfft excludes Nyquist point
+    S(2:end, :) = S(2:end, :).*2;
+else                                % even nfft includes Nyquist point
+    S(2:end-1, :) = S(2:end-1, :).*2;
+end
+
+% convert amplitude spectrum to dB (min = -120 dB)
+S = 20*log10(S + 1e-6);
+
+% plot the spectrogram
+figure(1)
+surf(t, f, S)
+shading interp
+axis tight
+box on
+view(0, 90)
+set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
+xlabel('Time, s')
+ylabel('Frequency, Hz')
+title('Amplitude spectrogram of the signal')
+
+handl = colorbar;
+set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
+ylabel(handl, 'Magnitude, dB')

Software/Matlab STFT Program/stft.m

@@ -0,0 +1,56 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%              Short-Time Fourier Transform            %
+%               with MATLAB Implementation             %
+%                                                      %
+% Author: M.Sc. Eng. Hristo Zhivomirov        12/21/13 %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [stft, f, t] = stft(x, wlen, hop, nfft, fs)
+
+% function: [stft, f, t] = stft(x, wlen, hop, nfft, fs)
+% x - signal in the time domain
+% wlen - length of the analysis Hamming window
+% hop - hop size
+% nfft - number of FFT points
+% fs - sampling frequency, Hz
+% stft - STFT matrix (only unique points, time across columns, freq across rows)
+% f - frequency vector, Hz
+% t - time vector, s
+
+% represent x as column-vector
+x = x(:);
+
+% length of the signal
+xlen = length(x);
+
+% form a periodic hamming window
+win = hamming(wlen, 'periodic');
+
+% stft matrix estimation and preallocation
+rown = ceil((1+nfft)/2);            % calculate the total number of rows
+coln = 1+fix((xlen-wlen)/hop);      % calculate the total number of columns
+stft = zeros(rown, coln);           % form the stft matrix
+
+% initialize the signal time segment index
+indx = 0;
+
+% perform STFT
+for col = 1:coln
+    % windowing
+    xw = x(indx+1:indx+wlen).*win;
+
+    % FFT
+    X = fft(xw, nfft);
+
+    % update the stft matrix
+    stft(:, col) = X(1:rown);
+
+    % update the index
+    indx = indx + hop;
+end
+
+% calculate the time and frequency vectors
+t = (wlen/2:hop:wlen/2+(coln-1)*hop)/fs;
+f = (0:rown-1)*fs/nfft;
+
+end