function [Y,MX] = logfsgram(X, N, SR, WIN, NOV, FMIN, BPO) % [Y,MX] = logfsgram(X, N, SR, WIN, NOV, FMIN, BPO) % Calculate a log-frequency spectrogram % X is input signal; N is parent FFT window; SR is the source samplerate. % WIN is actual window length within FFT, NOV is number of overlapping % points between successive windows. % Optional FMIN is the lowest frequency to display (80Hz); % BPO is the number of bins per octave (12). % MX returns the nlogbin x nfftbin mapping matrix; % sqrt(MX'*(Y.^2)) is an approximation to the original FFT % spectrogram that Y is based on, suitably blurred by going % through the log-F domain. % 2004-03-30 dpwe@ee.columbia.edu \$Header: /homes/dpwe/matlab/columbiafns/RCS/logfsgram.m,v 1.3 2004/04/01 22:39:02 dpwe Exp \$ if nargin < 2 N = 1024; end if nargin < 3 SR = 8000; end if nargin < 4 WIN = []; end if nargin < 5 NOV = []; end if nargin < 6 FMIN = 80; end if nargin < 7 BPO = 12; end if isempty(WIN) WIN = N; end if isempty(NOV) NOV = WIN/2; end % Calculate underlying STFT XX = specgram(X,N,SR,WIN,NOV); % Construct mapping matrix % Ratio between adjacent frequencies in log-f axis fratio = 2^(1/BPO); % How many bins in log-f axis nbins = floor( log((SR/2)/FMIN) / log(fratio) ); % Freqs corresponding to each bin in FFT fftfrqs = [0:(N/2)]*(SR/N); nfftbins = N/2+1; % Freqs corresponding to each bin in log F output logffrqs = FMIN * exp(log(2)*[0:(nbins-1)]/BPO); % Bandwidths of each bin in log F logfbws = logffrqs * (fratio - 1); % .. but bandwidth cannot be less than FFT binwidth logfbws = max(logfbws, SR/N); % Controls how much overlap there is between adjacent bands ovfctr = 0.5475; % Adjusted by hand to make sum(mx'*mx) close to 1.0 % Weighting matrix mapping energy in FFT bins to logF bins % is a set of Gaussian profiles depending on the difference in % frequencies, scaled by the bandwidth of that bin freqdiff = ( repmat(logffrqs',1,nfftbins) - repmat(fftfrqs,nbins,1) )./repmat(ovfctr*logfbws',1,nfftbins); mx = exp( -0.5*freqdiff.^2 ); % Normalize rows by sqrt(E), so multiplying by mx' gets approx orig spec back mx = mx ./ repmat(sqrt(2*sum(mx.^2,2)), 1, nfftbins); % Perform mapping in magnitude-squared (energy) domain y = sqrt( mx * (abs(XX).^2) ); % so, we lost phase information... if nargout < 1 imagesc([0 length(X)/SR],[1 nbins],20*log10(y)); axis xy xlabel('Time'); ylabel('Frequency'); yt = get(gca,'YTick'); for i = 1:length(yt) ytl{i} = sprintf('%.0f',logffrqs(yt(i))); end set(gca,'YTickLabel',ytl); else Y = y; MX = mx; end