spectrumAnalyzer function for DST (Matlab)

spectrumAnalyzer.m

% Michael Bartnett
% mrb402
% Friday November 4, 2011
% spectrumAnalyzer.m
%
% Function spectrumAnalyzer
%
% Reads a wav file and generates a STFT
%
% Usage: spectrumAnalyzer(filename, winLength, overlapLength, winType[, fftLength])
%
% filename - The wav file to read
%
% winLength - The length of each FFT window
%
% overlapLength - The amount to overlap between FFT windows
%
% winType - The windowing function to use
%
% fftLength - Optional parameter to increase the resolution of each FFT window
% 
% Window functions:
%
%   rect     - Uses rectwin()
%   hamming  - Uses hamming()
%   hann     - Uses hann()
%   blackman - Uses blackman()
%   bartlett - Uses bartlett
%   window   - Uses bartlett()
%
function y = spectrumAnalyzer(filename, winLength, overlapLength, winType, fftLength)
	if nargin < 4
		error('Parameters: filename, winLength, overlapLength, winType[, fftLength]');
	elseif nargin == 4
		fftLength = winLength;
	elseif (nargin == 5) && (fftLength <= winLength)
		error('If specified, fftLength must be greater than winLength')
	end

	if overlapLength > winLength
		error('winLength must be greater than overlapLength')
	end

	% Generate window based on winType parameter
	switch winType
	case 'rect'
		window = rectwin(winLength);
	case 'hamming'
		window = hamming(winLength);
	case 'hann'
		window = hann(winLength);
	case 'blackman'
		window = blackman(winLength);
	case 'bartlett'
		window = bartlett(winLength);
	otherwise
		error('winType must be one of rect, hamming, hann, blackman, or bartlett');
	end

	% Read in wav file, capture numchannels
	samples = wavread(filename);
	[numSamples, channels] = size(samples);

	% Sum to mono if more than one channel
	mix = samples(:,1);
	if channels > 1
		for i = 2:channels
			mix = mix + samples(:,i);
		end
	end

	% Normalize mix
	mix = mix / max(mix);

	% Calculate the number of frames, and how much to pad the source signal
	numFrames = ceil(numSamples / (winLength - overlapLength));
	targetLen = numFrames * (winLength - overlapLength) + overlapLength;;
	padCount = targetLen - length(mix);

	% Pad it
	mix = [mix; zeros(padCount,	 1)];

	% Preallocate STFT buffer
	y = zeros(numFrames, winLength);

	for i = 1:numFrames - 1
		r1 = i * (winLength - overlapLength) + 1; % calculate start of slice
		r2 = (r1 + winLength - 1);                % calculate end of slice

		% Pull slice out of source signal
		submix = mix(r1:(r1 + winLength - 1),1) .* window;


		if fftLength > winLength
			pre = ceil((fftLength - winLength) / 2);
			post = fftLength - pre;
			submix = [zeros(pre,1); submix; zeros(post,1)];
		end

		% Get the magnitude
		transformed = abs(fft(submix, winLength));

		% Add it to the next row in the preallocated STFT buffer
		y(i,:) = transformed;
	end

	% Normalize the magnitudes
	y = y / sum(window);

	% Scale magnitudes to dB
	y = log10(y / 1.0);
end