# FFT function isn't working on imported signals greater than a second

I have an fft function in matlab/octave that works great on an imported mono wav/audio signal that is 1 second long, but when the signal imported is more than 1 second the frequency values become incorrect. To check this I created a 2 second signal using audacity the first half of the 2 second signal (1 second of it) was 500hz at 0.3 amplitude and the second half of the 2 second signal (1 second of it) was 200hz at 0.8 amplitude. I also checked it with audacity's analyse spectrum option.

When I use the function on this 2 second signal it shows that the max frequency is at 999.95hz when it should show the max frequency is at 200hz. Can anyone help me get this function to work with signals greater than just 1 second.

I've included the function along with the code that calls the function.

The FUNCTION:

function [freq,amp,ampinv,phase,phase_radtodeg,phase_degtorad,phase_adjindeg,phase_adjinrad,sigfft,sigifft,sigphase]=rtfftphase(vp_sig_orig,Fs,phase_goal)

vp_sig_orig = vp_sig_orig - mean(vp_sig_orig); %remove the mean of your signal due to dc offset

vp_sig_orig=vp_sig_orig';
vp_sig_len=length(vp_sig_orig); %get sample rate from vp fs_rate needs to be an even number?

% Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
nfft= 2^(nextpow2(length(vp_sig_orig)));

% Take fft, padding with zeros so that length(fftx) is equal to nfft
fftx = fft(vp_sig_orig,nfft);
sigfft= fft(vp_sig_orig);
sigifft=ifft(sigfft);
sigphase = unwrap(angle(sigfft')); %get phase of orginal signal

% Calculate the number of unique points
NumUniquePts = ceil((nfft+1)/2);

% FFT is symmetric, throw away second half
fftx = fftx(1:NumUniquePts);

% Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
mx = abs(fftx)/length(vp_sig_orig); %replaced for testing from stackexchage

% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.
if rem(nfft, 2) % odd nfft excludes Nyquist point
mx(2:end) = mx(2:end)*2;
else
mx(2:end -1) = mx(2:end -1)*2;
end

%yamp=(mx(1,:)/max(abs(mx(1,:)))*1); % keep at 1, dont use in function use in script file directly
amp=mx;
ampinv=abs(amp-max(amp));

% This is an evenly spaced frequency vector with NumUniquePts points.
freq_vect = (0:NumUniquePts-1)*vp_sig_len/nfft;
freq=(freq_vect'); %take of round if you want exact numbers

%get phase of new signal
phase = unwrap(angle(fftx)); %get phase of orginal signal

%phase stuff phase adj stuff not fully done yet I don't know the correct way to use it leave at 0 for now


Example: that calls the function

imported_sig_L=imported_sig(:,1)';

%Sort Array Section

array1_sort_amp_norm=sortrows(array1,-2); %sort by amplitude
array1_sort_amp_inv=sortrows(array1,-3); %sort by amplitude inverse

%get max freq checking to make sure freq isn't zero if so go to next row
if (array1_sort_amp_norm(1,1)>0); % To get max amplitude
maxfreq_normA=array1_sort_amp_norm(1,1)
elseif   (array1_sort_amp_norm(1,1)<=0); %
maxfreq_normA=array1_sort_amp_norm(2,1)
end

if (array1_sort_amp_inv(1,1)>0); % To get inverted max amplitude
maxfreq_invA=array1_sort_amp_inv(1,1)
elseif   (array1_sort_amp_inv(1,1)<=0); %
maxfreq_invA=array1_sort_amp_inv(2,1)
end


Thanks

• I think that you forgot to divide frequency to sample length. Mar 23 '13 at 6:42
• @Eddy_Em I tried your suggestion and did freq=(freq_vect')./length(vp_sig_orig) unfortunately it didn't work instead of the freq being 200hz it was 0.0045319hz any other ideas? Mar 23 '13 at 7:00
• @Eddy_Em I thought the lines nfft= 2^(nextpow2(length(vp_sig_orig))); and freq_vect = (0:NumUniquePts-1)*vp_sig_len/nfft; would fix this I tried your suggestion and did freq=(freq_vect')./length(vp_sig_orig) unfortunately it didn't work instead of the freq being 200hz it was 0.0045319hz any other ideas? Mar 23 '13 at 7:46
• Look into my answer. Let's say your initial data have $N$ points and its length is $T$ seconds, then difference in neighbour points by T-axe would be $\Delta T = T / N$. In frequency domain you will have the same $N$ points, but the max frequency would be equal to $1 / \Delta T$ i.e. $N / T$. So, you just should divide $n$ (number of point of FFT) to $T$ (time length of data sample). Mar 23 '13 at 7:54

Here a short example of how to deal with FFT in octave/matlab:

function freq(npts, len, basefreq)
x = [0 : npts - 1] / npts * len;
y = sin(basefreq * x * 2 * pi);
F = abs(fft(y));
half = npts / 2;
F = F(1 : half);
plot([0 : half - 1] / len, F);
printf("Base frequency = %g\n", (find(F == max(F)) - 1) / len);
endfunction


try this function with different parameters.

Here npts is size of data vector; len - time length of data (in seconds); basefreq - base frequency for sine function. First we fit a time axe (x), after that we calculate sinusoid with given frequency, get modulus of its FFT, trim it to half (BTW, another variant is not to trim, but to move frequency zero to middle of fftshift'ed data). After all we plot result of FFT and print out founded base frequency.

• Thanks this helped a lot. I just need to alter one line freq=(freq_vect').*Fs/vp_sig_len; Mar 23 '13 at 22:20