Dawn@4: function [yo,fo,to] = h_specgram2(varargin) Dawn@4: %SPECGRAM Calculate spectrogram from signal. Dawn@4: % B = SPECGRAM(A,NFFT,Fs,WINDOW,SHIFT) calculates the spectrogram for Dawn@4: % the signal in vector A. SPECGRAM splits the signal into overlapping Dawn@4: % segments, windows each with the WINDOW vector and forms the columns of Dawn@4: % B with their zero-padded, length NFFT discrete Fourier transforms. Thus Dawn@4: % each column of B contains an estimate of the short-term, time-localized Dawn@4: % frequency content of the signal A. Time increases linearly across the Dawn@4: % columns of B, from left to right. Frequency increases linearly down Dawn@4: % the rows, starting at 0. If A is a length NX complex signal, B is a Dawn@4: % complex matrix with NFFT rows and Dawn@4: % k = fix((NX-NOVERLAP)/(length(WINDOW)-NOVERLAP)) Dawn@4: % columns, where NOVERLAP = length(WINDOW)-mean(SHIFT) Dawn@4: % If A is real, B still has k columns but the higher frequency Dawn@4: % components are truncated (because they are redundant); in that case, Dawn@4: % SPECGRAM returns B with NFFT/2+1 rows for NFFT even and (NFFT+1)/2 rows Dawn@4: % for NFFT odd. If you specify a scalar for WINDOW, SPECGRAM uses a Dawn@4: % Hanning window of that length. WINDOW must have length smaller than Dawn@4: % or equal to NFFT and greater than NOVERLAP. NOVERLAP is the number of Dawn@4: % samples the sections of A overlap. Fs is the sampling frequency Dawn@4: % which does not effect the spectrogram but is used for scaling plots. Dawn@4: % Dawn@4: % [B,F,T] = SPECGRAM(A,NFFT,Fs,WINDOW,NOVERLAP) returns a column of Dawn@4: % frequencies F and one of times T at which the spectrogram is computed. Dawn@4: % F has length equal to the number of rows of B, T has length k. If you Dawn@4: % leave Fs unspecified, SPECGRAM assumes a default of 2 Hz. Dawn@4: % Dawn@4: % B = SPECGRAM(A) produces the spectrogram of the signal A using default Dawn@4: % settings; the defaults are NFFT = minimum of 256 and the length of A, a Dawn@4: % Hanning window of length NFFT, and NOVERLAP = length(WINDOW)/2. You Dawn@4: % can tell SPECGRAM to use the default for any parameter by leaving it Dawn@4: % off or using [] for that parameter, e.g. SPECGRAM(A,[],1000) Dawn@4: % Dawn@4: % See also PWELCH, CSD, COHERE and TFE. Dawn@4: Dawn@4: % Author(s): L. Shure, 1-1-91 Dawn@4: % T. Krauss, 4-2-93, updated Dawn@4: % Copyright 1988-2000 The MathWorks, Inc. Dawn@4: % $Revision: 1.6 $ $Date: 2000/06/09 22:07:35 $ Dawn@4: Dawn@4: error(nargchk(1,5,nargin)) Dawn@4: [msg,x,nfft,Fs,window,shift]=specgramchk(varargin); Dawn@4: error(msg) Dawn@4: Dawn@4: nx = length(x); Dawn@4: nwind = length(window); Dawn@4: if nx < nwind % zero-pad x if it has length less than the window length Dawn@4: x(nwind)=0; nx=nwind; Dawn@4: end Dawn@4: x = x(:); % make a column vector for ease later Dawn@4: window = window(:); % be consistent with data set Dawn@4: noverlap = nwind - mean(shift); Dawn@4: ncol = fix((nx-noverlap)/(nwind-noverlap)); Dawn@4: Dawn@4: totalshift = sum(shift); Dawn@4: numshift = length(shift); Dawn@4: progshift = (cumsum(shift)-shift(1))'; Dawn@4: overshift = 1:totalshift:nx; Dawn@4: shifts = progshift(:,ones(1,length(overshift)))+overshift(ones(length(progshift),1),:); Dawn@4: colindex = reshape(shifts,1,prod(size(shifts))); Dawn@4: colindex = colindex(1:ncol); Dawn@4: %colindex = 1 + (0:(ncol-1))*(nwind-noverlap); Dawn@4: rowindex = (1:nwind)'; Dawn@4: if length(x)<(nwind+colindex(ncol)-1) Dawn@4: x(nwind+colindex(ncol)-1) = 0; % zero-pad x Dawn@4: end Dawn@4: Dawn@4: y = zeros(nwind,ncol); Dawn@4: Dawn@4: % put x into columns of y with the proper offset Dawn@4: % should be able to do this with fancy indexing! Dawn@4: y(:) = x(rowindex(:,ones(1,ncol))+colindex(ones(nwind,1),:)-1); Dawn@4: Dawn@4: % Apply the window to the array of offset signal segments. Dawn@4: y = window(:,ones(1,ncol)).*y; Dawn@4: Dawn@4: % now fft y which does the columns Dawn@4: y = fft(y,nfft); Dawn@4: if ~any(any(imag(x))) % x purely real Dawn@4: if rem(nfft,2), % nfft odd Dawn@4: select = [1:(nfft+1)/2]; Dawn@4: else Dawn@4: select = [1:nfft/2+1]; Dawn@4: end Dawn@4: y = y(select,:); Dawn@4: else Dawn@4: select = 1:nfft; Dawn@4: end Dawn@4: f = (select - 1)'*Fs/nfft; Dawn@4: Dawn@4: t = (colindex-1)'/Fs; Dawn@4: Dawn@4: % take abs, and use image to display results Dawn@4: if nargout == 0 Dawn@4: % newplot; Dawn@4: % if length(t)==1 Dawn@4: % imagesc([0 1/f(2)],f,20*log10(abs(y)+eps));axis xy; colormap(jet) Dawn@4: % else Dawn@4: % imagesc(t,f,20*log10(abs(y)+eps));axis xy; colormap(jet) Dawn@4: % end Dawn@4: % xlabel('Time') Dawn@4: % ylabel('Frequency') Dawn@4: elseif nargout == 1, Dawn@4: yo = y; Dawn@4: elseif nargout == 2, Dawn@4: yo = y; Dawn@4: fo = f; Dawn@4: elseif nargout == 3, Dawn@4: yo = y; Dawn@4: fo = f; Dawn@4: to = t; Dawn@4: end Dawn@4: Dawn@4: function [msg,x,nfft,Fs,window,shift] = specgramchk(P) Dawn@4: %SPECGRAMCHK Helper function for SPECGRAM. Dawn@4: % SPECGRAMCHK(P) takes the cell array P and uses each cell as Dawn@4: % an input argument. Assumes P has between 1 and 5 elements. Dawn@4: Dawn@4: msg = []; Dawn@4: Dawn@4: x = P{1}; Dawn@4: if (length(P) > 1) & ~isempty(P{2}) Dawn@4: nfft = P{2}; Dawn@4: else Dawn@4: nfft = min(length(x),256); Dawn@4: end Dawn@4: if (length(P) > 2) & ~isempty(P{3}) Dawn@4: Fs = P{3}; Dawn@4: else Dawn@4: Fs = 2; Dawn@4: end Dawn@4: if length(P) > 3 & ~isempty(P{4}) Dawn@4: window = P{4}; Dawn@4: else Dawn@4: if length(nfft) == 1 Dawn@4: window = hanning(nfft); Dawn@4: else Dawn@4: msg = 'You must specify a window function.'; Dawn@4: end Dawn@4: end Dawn@4: if length(window) == 1, window = hanning(window); end Dawn@4: if (length(P) > 4) & ~isempty(P{5}) Dawn@4: shift = P{5}; Dawn@4: else Dawn@4: shift = ceil(length(window)/2); Dawn@4: end Dawn@4: Dawn@4: % NOW do error checking Dawn@4: if (length(nfft)==1) & (nfft