1
$\begingroup$

I am new to matlab. I have a noisy dataset represented as a signal on which I have to compute the fft and plot the one-sided spectrum, however I'm quite unsure if my data's signal length I assume that is denoted by the size of my vector (4x10000double).

Just a picture of a piece of how my data looks:

Data

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal  
t = (0:L-1)*T;        % Time vector

This ^ is what I see on almost every implementation, and I do not understand if there is a way to adapt it to my specific case (for example L=40000 in my case?) and then I clearly do not understand how to compute my sampling frequency, or how would I change my approach of computing the spectra.

In the guide I also see the following computations for the double, and then single-sided spectra of the signal after fft is applied:

Y = fft(X); // compute the fft

P2 = abs(Y/L); //- compute double-sided spectra by abs() of the DFT values over length of signal (unsure about how this works)

P1 = P2(1:L/2+1); //here is the single sided! this one I get...

P1(2:end-1) = 2*P1(2:end-1); //what is this?

f = Fs*(0:(L/2))/L; //a bit unsure about this as well, shouldn't it be 0:500 all the time?

plot(f,P1) // just ploting

Help on this would be more than apreciated! I sincerely need it :D

$\endgroup$

2 Answers 2

1
$\begingroup$

You can do the Fourier transform without knowing the sampling frequency, in the same way you could plot that data, the only difference is that the X axis is in arbitrary units, instead of something meaningful.

It seems clear to me that you are not yet very familiar with the Fourier transform (that is completely OK). However, you should first understand why you want to perform a Fourier transform and whether you do need to know the sampling frequency for your purposes or not. There are cases where you don't really care about it, and sometimes doing a Fourier transform is just an intermediate step in an operation (where you are going to do the inverse Fourier transform after).

Now to teach you some practical tidbits about the FFT:

Y = fft(X); // compute the fft P2 = abs(Y/L); //- compute double-sided spectra by abs() of the DFT values over length of signal (unsure about how this works)

when you do the fft (the fast Fourier transform algorithm) you get the full, complex, Fourier transform scaled by the length of your vector (L). You divide it by L to get the correct amplitude.

The reason you do abs(...) is because you want to plot the amplitude of the frequency components that make up your signal, but not their phase. The discrete Fourier transform gives you a vector of complex values (they have a real and imaginary part, or in a more instructive way to put it, a magnitude and a phase). For real signals, each corresponds to a sine wave with the amplitude given by the abs(...) and the phase given by the angle(...).

P1 = P2(1:L/2+1); //here is the single sided! this one I get... P1(2:end-1) = 2*P1(2:end-1); //what is this?

When you do the FFT for real signals, the left and right side contain the same information, but half of the amplitude in in each side. So, after getting only one of the single sides, you multiply everything by 2. The only exception is the first value of the vector (the zero frequency, DC offset, bias, or whatever you prefer calling it). That component is common to both of the "sides" of the FFT, so you don't double it.

f = Fs*(0:(L/2))/L; //a bit unsure about this as well, shouldn't it be 0:500 all the time?

This is to generate the correct X vector for your plot (i.e., the frequency vector). This is what you need the sampling frequency for. The length of the X vector after the Fourier transform (frequency vector) is the same as the length of the X vector before the Fourier transform (time vector). Since you are doing single sided, half of the values go away.

Essentially what that line of code does is get you a vector with L/2 values that goes from 0 to 0.5. Then you multiply that by the sampling frequency and you get a vector of equally spaced values that go from 0 to Fs/2, or the Nyquist frequency.

plot(f,P1) // just ploting

$\endgroup$
1
  • $\begingroup$ Nice answer! +1 from me. I was in the middle of writing my answer when yours posted so I figured I might as well keep going. Both say the same thing with different wording, which will be beneficial for the OP. $\endgroup$
    – Jdip
    Feb 28, 2023 at 22:48
1
$\begingroup$

Ok let's take this one step at a time:

  1. Let's call $x(t)$ the continuous signal that was sampled at some unknown sampling frequency $f_s = \frac{1}{T_s}$ to give you the $N = \mathtt{40000}$ samples $x[n] = x(nT_s)$ in your dataset.

  2. The Discrete Fourier Transform ($\texttt{DFT}$) of $x[n]$ is defined as: $$X[k] = \sum_{n=0}^{N-1}x[n]e^{-i\frac{2\pi}{N}kn} \qquad k = 0\dots N-1 \tag{1}$$ Notice there is no mention of sampling frequency $f_s$ here. The DFT is a transform, it does not care about sampling frequency. Where the sampling frequency comes in play is a little later down this chain of steps.
    One thing to know is for real input $x[n]$, $X[k]$ is complex conjugate symmetric: $$X[k] = X^*[-k]$$ In plain English, this means that half of $X[k]$ is symmetric to the other half (even symmetry for the magnitude $|X[k]|$ and odd symmetry for the angle $\angle X[k]$).

  3. The $\texttt{DFT}$ is computed through the Fast Fourier Transform ($\texttt{FFT}$) algorithm: X = fft(x)/N
    Why do we scale by $N$? Here's one way to look at it: $X[0]$ gives you the fourier coefficient for the $0$-frequency, meaning the $\texttt{DC}$ component of your signal, i.e the average. If you compute it with $(1)$, you get $X[0] = \sum_{n=0}^{N-1}x[n]$. To get the average, you'd want to divide by $N$ right? Same for all other fourier coefficients ($k = 1,2,\dots,N-1$)

  4. Take the magnitude (X was scaled at the previous step): abs(X)

  5. Take one side since abs(X) is even-symmetric and one half is redundant: X = X(1:N/2+1)

  6. The $\texttt{DFT}$ splits the amplitude of each coefficient between the $2$ symmetric parts, so multiply by $2$, except for the $\texttt{DC}$ component (since $X[0]$ and $X^*[0]$ combine): X = 2*X(2:end-1)

  7. As promised, here is where the sampling frequency $f_s$ comes in play. At this point, each $\texttt{DFT}$ coefficient $X[k]$ is a value corresponding to a bin number $k = 0,\dots,N-1$. We need to match these bin numbers to a physical quantity (frequency in $\texttt{Hz}$) to plot $X[k]$ on the y-axis, and frequency on the x-axis. That's where we use $f_s$ to map each bin number $k$ to a frequency $f_k$ in $\texttt{Hz}$:
    $$f_k = k\frac{f_s}{N}$$ In code, this is done with f = fs/N*(0:(N/2));

  8. Plot! plot(f, X)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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