annotate utilities/UTIL_FFT.m @ 38:c2204b18f4a2 tip

End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents
children
rev   line source
rmeddis@38 1 function [fft_powerdB, fft_phase, frequencies, fft_ampdB]=...
rmeddis@38 2 UTIL_FFT(getFFTsignal, dt, referenceAmplitude)
rmeddis@38 3 % UTIL_FFT
rmeddis@38 4 % [fft_powerdB, fft_phase, frequencies, fft_ampdB]= UTIL_FFT(getFFTsignal, dt, referenceAmplitude)
rmeddis@38 5
rmeddis@38 6 x=length(getFFTsignal); % adjust the length of the signal to a power of 2 for FFT
rmeddis@38 7 n=2;
rmeddis@38 8 while n<=x
rmeddis@38 9 n=n*2;
rmeddis@38 10 end
rmeddis@38 11 n=round(n/2);
rmeddis@38 12 getFFTsignal=getFFTsignal(1:n);
rmeddis@38 13 frequencies = (1/dt)*(1:n/2)/n;
rmeddis@38 14
rmeddis@38 15 fft_result = fft(getFFTsignal, n); % Compute FFT of the input signal.
rmeddis@38 16 fft_power = fft_result .* conj(fft_result); % / n; % Compute power spectrum.
rmeddis@38 17 % Dividing by 'n' we get the power spectral density.
rmeddis@38 18 % fft_power = fft_result .* conj(fft_result) / n; % Compute power spectralDensity.
rmeddis@38 19
rmeddis@38 20 % Prepare the FFT for display
rmeddis@38 21 fft_power=fft_power(1:length(fft_power)/2); % remove mirror frequencies
rmeddis@38 22 % fft_powerdB = UTIL_amp2dB (fft_power, max(fft_power)); % convert to dB
rmeddis@38 23 fft_powerdB = 10*log10(fft_power); % convert to dB
rmeddis@38 24
rmeddis@38 25 % amplitude spectrum
rmeddis@38 26 if nargin<3
rmeddis@38 27 referenceAmplitude=28e-6; % for SPL
rmeddis@38 28 end
rmeddis@38 29 % refAmpdB= referenceAmplitude^0.5;
rmeddis@38 30 fft_ampdB = 10*log10(fft_power/referenceAmplitude); % convert to dB
rmeddis@38 31 % peak_fft_powerdBSPL=max(fft_ampdB)
rmeddis@38 32
rmeddis@38 33
rmeddis@38 34 fft_phase = angle(fft_result); % Compute the phase spectrum.
rmeddis@38 35 fft_phase=fft_phase(1:length(fft_phase)/2); % remove mirror phases
rmeddis@38 36 jags=find(diff(fft_phase)>0); % unwrap phase
rmeddis@38 37 for i=jags,
rmeddis@38 38 fft_phase(i+1:end)=fft_phase(i+1:end)-2*pi;
rmeddis@38 39 end