0
$\begingroup$

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
    phase_radtodeg=rtrad2deg(phase); %also check sigphase to see if it's the same.
    phase_degtorad=rtdeg2rad(phase_radtodeg); %for checking
    phase_adjindeg=phase_goal-phase_radtodeg; %used to create new phase to get to goal phase
    phase_adjinrad=rtdeg2rad(phase_adjindeg);

Example: that calls the function

imported_sig_L=imported_sig(:,1)';
[xfreq_orig,yamp_orig,yamp_inv,phase_orig,phase_in_deg,phase_in_rad,phase_adjindeg,phase_adjinrad,sigfft,sigifft,sigphase]=rtfftphase(imported_sig_L,fs_rate,0); 


array1=[xfreq_orig,yamp_orig,yamp_inv,phase_orig,mod(phase_in_deg,360),phase_in_rad,phase_adjindeg,phase_adjinrad]; %have degs only go up to 360


%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

$\endgroup$
  • $\begingroup$ I think that you forgot to divide frequency to sample length. $\endgroup$ – Eddy_Em Mar 23 '13 at 6:42
  • $\begingroup$ @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? $\endgroup$ – Rick T Mar 23 '13 at 7:00
  • $\begingroup$ @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? $\endgroup$ – Rick T Mar 23 '13 at 7:46
  • 1
    $\begingroup$ 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). $\endgroup$ – Eddy_Em Mar 23 '13 at 7:54
2
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.