rmeddis@38: function [fft_powerdB, fft_phase, frequencies, fft_ampdB]=...
rmeddis@38:     UTIL_FFT(getFFTsignal, dt, referenceAmplitude)
rmeddis@38: % UTIL_FFT
rmeddis@38: % [fft_powerdB, fft_phase, frequencies, fft_ampdB]= UTIL_FFT(getFFTsignal, dt, referenceAmplitude)
rmeddis@38: 
rmeddis@38: x=length(getFFTsignal);         % adjust the length of the signal to a power of 2 for FFT
rmeddis@38: n=2;
rmeddis@38: while n<=x
rmeddis@38:     n=n*2;
rmeddis@38: end
rmeddis@38: n=round(n/2);
rmeddis@38: getFFTsignal=getFFTsignal(1:n);
rmeddis@38: frequencies = (1/dt)*(1:n/2)/n;
rmeddis@38: 
rmeddis@38: fft_result = fft(getFFTsignal, n);				    % Compute FFT of the input signal.
rmeddis@38: fft_power = fft_result .* conj(fft_result); % / n;	% Compute power spectrum.  
rmeddis@38: % Dividing by 'n' we get the power spectral density.
rmeddis@38: %     fft_power = fft_result .* conj(fft_result) / n;	% Compute power spectralDensity.  
rmeddis@38: 
rmeddis@38: % Prepare the FFT for display
rmeddis@38: fft_power=fft_power(1:length(fft_power)/2);         % remove mirror frequencies
rmeddis@38: %     fft_powerdB = UTIL_amp2dB (fft_power, max(fft_power)); % convert to dB
rmeddis@38: fft_powerdB = 10*log10(fft_power);                  % convert to dB
rmeddis@38: 
rmeddis@38: % amplitude spectrum
rmeddis@38: if nargin<3
rmeddis@38:     referenceAmplitude=28e-6; % for SPL
rmeddis@38: end
rmeddis@38: % refAmpdB= referenceAmplitude^0.5;
rmeddis@38: fft_ampdB = 10*log10(fft_power/referenceAmplitude); % convert to dB
rmeddis@38: %     peak_fft_powerdBSPL=max(fft_ampdB)
rmeddis@38: 
rmeddis@38: 
rmeddis@38: fft_phase = angle(fft_result);			            % Compute the phase spectrum.
rmeddis@38: fft_phase=fft_phase(1:length(fft_phase)/2);         % remove mirror phases
rmeddis@38: jags=find(diff(fft_phase)>0);                       % unwrap phase
rmeddis@38: for i=jags, 
rmeddis@38:     fft_phase(i+1:end)=fft_phase(i+1:end)-2*pi; 
rmeddis@38: end