Fourier transform of $\cos(n\omega t)$

Question

My question is probably very stupid, but I've been strugling for a while on it now...

In need to find the Fourier transform of $1+\cos^3(2\pi ft)$.

I wrote that : $$\cos^3(2\pi ft)=\frac{\cos(6\pi ft)+3\cos(2\pi ft)}{4}$$ And so I have: $$\delta(f) +\frac 18 \bigg[\delta(f-3f_0)+\delta(f-3f_0)\bigg] +\frac 38 \bigg[\delta(f-f_0)+\delta(f-f_0)\bigg]$$

So, on my spectrum, I should have a dirac at $0$, a smaller one at $f_0$ and a smaller at $3f_0$...

But when I process it with matlab (using fast fourier transform), I get this :

(With a frequency of $10\textrm{ kHz}$).

So the dirac I thought would be at $3f_0$ is in fact at $\frac{f_0}{2}$. What am I missing ?

Answer 1

Let's take a look at the first half of your expansion; $cos(6\pi f_{0} t)$

The Fourier transform for this would be

\begin{equation} X_{c}(j\Omega) = \pi \delta(\Omega - 6\pi f_{0}) + \pi \delta(\Omega + 6\pi f_{0}) \end{equation}

For your Fourier transform to be correct, we need that \begin{equation} 6\pi f_{0} < \pi f_{s} \end{equation}

This means that your sampling rate must be high enough to avoid aliasing. Are you sure that your $f_{0} < f_{s}/2?$

To me it looks like this is the problem.

Answer 2

Thank you so much everyone!

I am really not confident in my mathematics skills and so I was focusing on the Fourier transform, but as you all guessed, the problem was my sampling frequency way too small (I was using 25kHz). With 60kHz it works like a charm. Still feels like an idiot though ;-)

Thank you !

PS : To sum up if anyone was facing the same kind of problem

The theoric part was wright. We have : $$\delta(f) +\frac 18 \bigg[\delta(f-3f_0)+\delta(f-3f_0)\bigg] +\frac 38 \bigg[\delta(f-f_0)+\delta(f-f_0)\bigg]$$

The problem was with my simulation: I forgot to check if my sampling frequency was meeting the Nyquist-Shannon criteria ($Fs > 2F_{max}$).

So, with this code :

Fs = 7e4;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 100000;             % Length of signal
t = (0:L-1)*T;        % Time vector
fm = 10e3;

X = 1+(cos(2*pi*fm*t)).^3;
%X = 1 + 0.25*cos(3*2*pi*10000*t)+(3/4)*cos(2*pi*10000*t);
Y = abs(fft(X));

f = 0:Fs/L:Fs/2;
plot(f,Y(1:L/2+1))
xlabel('f (Hz)')
ylabel('|P1(f)|')

We get the expected result: