% FP: function to compute & plot amplitude spectrum
% Usage: [a,f]=fp(wave,rate,window);
% wave: input waveform
% rate: sample rate in Hz (default 10000 Hz)
% window options: 'hann', 'hamm', 'kais', or 'rect' (default=hamming)
% [a,f]: log magnitude (dB re:1), frequency

function [a,f]=fp(x,rate,window);

if ~exist('rate','var'), 
    rate=10000; 
end;
if ~exist('window','var'),
    window='hamm'; 
end;

x=x(:);  % make x a column vector
n=length(x);

% apply window function
window=lower(window);
if window=='rect',
    x=x.*ones(n,1);
elseif window=='hamm',
    x=x.*hamming(n);
elseif window=='hann',
    x=x.*hanning(n);
elseif window=='kais',
    x=x.*kwind(ones(n,1),n/2-1);
end;
% compute log magnitude spectrum
m=fft(x,n);
no2=round(n/2);
a=20*log10( abs ( m ) / n);  
f=rate/n*(0:no2)/1000;
freq=f(1:no2);
amp=a(1:no2);
% make plot
plot(freq,amp);
xlabel('Frequency (kHz)');
ylabel('Amplitude (dB)');
axis([0 rate/2000 -Inf Inf]);