3.11 Plot the power spectral density of a signal

Problem: Given a signal \[ 3\sin \left ( 2\pi f_{1}t\right ) +4\cos \left ( 2\pi f_{2}t\right ) +5\cos \left ( 2\pi f_{3}t\right ) \] for the following frequency values (in Hz) \[ f_{1}=20 ,f_{2}=30, f_{3}=50 \] and plot the signal power spectral density so that 3 spikes will show at the above 3 frequencies. Use large enough sampling frequency to obtain a sharper spectrum

3.11.1 Mathematica

Clear ["Global`*"]; 
f1 = 20; 
f2 = 30; 
f3 = 100; 
maxTime = 0.2; 
 
y[t_]:=2 Sin[2 Pi f1 t]+4 Cos[2 Pi f2 t]+ 
       5 Cos[2 Pi f3 t]; 
 
Plot[y[t],{t,0,maxTime}, 
  Frame->True, 
  FrameLabel->{{"y(t)",None}, 
      {"t (sec)","signal in time domain"}}, 
  RotateLabel->False, 
  GridLines->Automatic, 
  GridLinesStyle->Dashed, 
  PlotRange->All,BaseStyle -> 12]
 

fs  = 7.0*Max[{f1,f2,f3}]; 
yn  = Table[y[n],{n,0,maxTime,(1/fs)}]; 
len = Length[yn]; 
py  = Fourier[yn]; 
nUniquePts = Ceiling[(len+1)/2]; 
py  = py[[1;;nUniquePts]]; 
py  = Abs[py]; 
py  = 2*(py^2); 
py[[1]] = py[[1]]/2; 
f   = (Range[0,nUniquePts-1] fs)/len; 
 
ListPlot[Transpose[{f,py}],Joined->False, 
  FrameLabel->{{"|H(f)|",None}, 
               {"hz","Magnitude spectrum"}}, 
  ImageSize->400, 
  Filling->Axis, 
  FillingStyle->{Thick,Red}, 
  Frame->True, 
  RotateLabel->False, 
  GridLines->Automatic, 
  GridLinesStyle->Dashed, 
  PlotRange->All,BaseStyle -> 12]
 

Here is the same plot as above, but using InterpolationOrder -> 0

And using InterpolationOrder -> 2

3.11.2 Matlab

clear all; close all; 
maxTime = 0.2; 
t  = 0:0.001:maxTime; 
f1 =20; %Hz 
f2 =30; %Hz 
f3 =100; %Hz 
y = @(t) 2*sin(2*pi*f1.*t)+4*cos(2*pi*f2.*t)+... 
         5*cos(2*pi*f3.*t); 
plot(t,y(t)); 
set(gcf,'Position',[10,10,320,320]); 
grid on 
title('signal in time domain'); 
xlabel('time(msec)'); 
ylabel('Amplitude');
 

fs  = 7.0*max([f1 f2 f3]); 
delT=1/fs; 
yn  = y(linspace(0,maxTime,maxTime/delT)); 
len = length(yn); 
py  = fft(yn); 
nUniquePts = ceil((len+1)/2); 
py  = py(1:nUniquePts); 
py  = abs(py); 
py  = 2*(py.^2); 
py(1) = py(1)/2; 
f   = 0:nUniquePts-1; 
f   = f*fs/len; 
plot(f,py); 
title('power spectral density'); 
xlabel('freq (Hz)'); 
ylabel('Energy content'); 
grid on 
set(gcf,'Position',[10,10,320,320]);