Mercurial > hg > map
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 |