function [Ba2,F,T,Bl10]=spectrogram(x,nfft,Fs,wlen,olap,unt) % [Ba2,F,T,Bl10]=spectrogram(x,nfft,Fs,wlen,olap,unt) % % Computation of TIME-DEPENDENT SPECTRAL DENSITY (UNIT^2/HZ) % by overlapping segment analysis with a Hanning window. % % INPUT: % % x Real-valued signal to be analyzed % nfft Number of frequencies calculated (default: 256) % Fs Sampling frequency (Hz) used only for output scaling (default: 1) % wlen Length of windowed segment, in samples, which must be % smaller than or equal to 'nfft' and larger than 'olap' % (default: 256) % olap Overlap of data segments, in samples (default: 70) % unt String with the unit name [default: s] % % OUTPUT: % % Ba2 Each column of Ba2 contains an estimate of the short-term, % time-localized spectral density of the signal. (UNITS^2/(1/unt)) % Time increases linearly across the columns of Ba2, from left to % right. Frequency increases linearly down the rows, starting % at 0. Ba2 is a real matrix with k columns, where % k = floor((length(x)-olap)/(wlen-olap)), and a number of rows % equal to nfft/2+1 if 'x' is real and 'nfft' even, % and (nfft+1)/2 for 'x' is real and 'nfft' odd. % Note that our Ba2 is abs(fft(x).^2)/Fs. % F Frequency axis (1/unt, which is Hz by default) % T Time axis (unt), starting from zero % Bl10 10*log10(Ba2) % % See also PCHAVE, TIMSPECPLOT, SPECTROGRAM2 % % EXAMPLE: % % Plot the results with, e.g.: % % IMAGESC(h1.B+wlen/Fs/2+T,F,Bl10); where h1.B is the begin time % % Last modified by fjsimons-at-alum.mit.edu, 10/22/2012 defval('wlen',256) defval('olap',70) defval('nfft',256) defval('Fs',1) defval('unt','s') if nfft>wlen disp(sprintf('SPECTROGRAM: NFFT of %i is larger than window length of %i',... nfft,wlen)) end if nfft=wlen, error('Overlap must be strictly smaller than window length') end disp(sprintf('Window size for SPECTROGRAM: %g %s',wlen/Fs,unt)) % OVERLAPPING SEGMENT ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compute Hanning window dwin=fhanning(wlen); % Normalize window so the sum of squares is unity (PW 208a) dwin=dwin/sqrt(dwin'*dwin); % Determine window parameters and truncation of data set npts=length(x); eflen=wlen-olap; checkit=(npts-olap)/eflen; nwin=floor(checkit); disp(sprintf('Number of overlapping data segments: %i',nwin)) if nwin~=checkit disp(sprintf(... 'Number of segments is not integer: %i / %i points truncated',... npts-(nwin*eflen+olap),npts)) end if nwin<=0; error('Data sequence not long enough'); end % Make matrix out of suitably repeated windowed segments % of the data 'xsdw' is x segmented THEN detrended THEN % windowed with normalized window rows=[1:wlen]; cols=[0:(nwin-1)]*eflen; xsd=detrend(... x(repmat(rows(:),1,nwin)+... repmat(cols,wlen,1))); xsdw=xsd.*repmat(dwin,1,nwin); % This in case it's a boxcar window! % xsdb=xsd/sqrt(wlen); % Calculate POWER SPECTRAL DENSITY (power per frequency) Ba2=abs(fft(xsdw,nfft,1)).^2; % Next thing is the SPECTRAL DENSITY (energy per frequency) thus % in UNITS^2/HZ or UNITS^2*SECOND Ba2=Ba2/Fs; % Collect frequency information and select output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculate frequency vector for REAL signals % and get rid of periodicity in spectrum selekt=[1:floor(nfft/2)+1]; F=(selekt-1)'*Fs/nfft; % Calculate time axis; cols is the sample number T=cols(:)/Fs; % Need to scale frequencies except 0 and Nyquist by a factor of two % if you only take half of the spectrum (for real signals) % Compute one-sided spectrum (Bendat and Piersol Eq. 5.33 and page 424.) Ba2=Ba2(selekt,:); Ba2=[Ba2(1,:); 2*Ba2(2:end-1,:); Ba2(end,:)]; Bl10=10*log10(Ba2);