Amplitude and phase spectrum in MATLAB

Question

I have a little problem with the subject Theory of signal. We were told to make amplitude and phase spectrum of cosinus in MATLAB, but we don't have results nor materials from our teachers, so we have to study from many different materials and you can imagine what this is like. So I'm a little confused.

However, I did some graphs, that should be what I want, according to instructions, but I have never seen something like that, so I'm not sure, if that's right.

Do these graphs make any sense to you? I used fast Fourier transform (fft(y) in MATLAB) on the signal and then plot the result of transform versus sample frequency and I thought it's possible to get phase of the signal by angle(fft(y)), so I plotted this on the last graph. I'm very unsure especially by the third graph.

I would like to apologize for incorrect notions I might use, 'cause I'm not a native english speaker.

Thanks for any help.

Answer 1

This is the most simple FT that you should understand. The basic rule is: a sine of frequency f0 in time domain is a delta at f0 and -f0 in frequency. You can do this theoretically, as shown in http://mathworld.wolfram.com/FourierTransformCosine.html.

Those peaks in your graph are consistent with this. However, it isn't true that there is a peak in at 90Hz. The peaks should be at 10 Hz and -10 Hz. There are two possible approaches here:

a) Use fftshift to correct this issue. This will take each frequency to its real place.

Fs=100; %sampling frequency
T=1; %signal length
N=T*Fs;     %number of samples
f=-Fs/2:Fs/N:Fs/2-Fs/N; %frequency vector
x=cos(2*pi*10*t);
Xf=fft(x);
Xf=fftshift(Xf);
plot(f,abs(Xf)) %magnitude

b) Assume the signal is real, so the spectrum is symmetric. It is true in this case, but it doesn't have to. You could just plot the first half of your fft, i.e. from 0 Hz to 50 Hz.

By the way, if you plot plot(f,abs(Xf)/N), you will find that the peaks have an amplitude of 0.5A (being A the amplitude of the sinusoid). The phase plot seems OK to me, especially after fftshift-ing it.

Answer 2

Look at this MATLAB function, it can calculate phase spectrum as well as amplitude spectrum with a perfect accuracy.

This program calculates amplitude and phase spectra of an input signal with acceptable accuracy especially in the calculation of phase spectrum.The code does three main jobs for calculation amplitude and phase spectra. First of all, it extends the input signal to infinity; because for calculation Fourier transform(FT) (fft function in Matlab), we consider our signal is periodic with an infinite wavelength, the code creates a super_signal by putting original signal next to itself until the length of super_signal is around 1000000 samples, why did I choose 1000000 samples? Actually, it is just based on try and error!! For me, a supper signal with 1000000 samples has the best output.

Second, for calculating fft in Matlab you can choose different resolutions, the Mathwork document and help use NFFT=2^nextpow2(length(signal)), it definitely isn't enough for one that wants high accuracy output. Here, I choose the resolution of NFFT=100000 that works for most signals.

Third, programs filters result of FT by thresholding, it is very important! For calculating phase spectrum, its result is very noisy because of floating rounding off error, it causes during calculation arctan even small rounding off error produces significant noise in the result of phase spectrum, for suppressing this kind of noise you can define a threshold value. It means if amplitude of specific frequency is less than a predefined threshold value (you must define it) it put zero instead of it.

These three steps help to improve the result of amplitude and phase spectra significantly.

IF YOU USE THIS PROGRAM IN YOUR RESEARCH, PLEASE CITE THE FOLLOWING PAPER:

Afshin Aghayan, Priyank Jaiswal, and Hamid Reza Siahkoohi (2016). ”Seismic denoising using the redundant lifting scheme.” GEOPHYSICS, 81(3), V249-V260. https://doi.org/10.1190/geo2015-0601.1

Thank you! Afshin Aghayan 2017

% This program calculates amplitude and phase spectrums of an input signal
% with acceptable accuracy especially in the calculation of phase spectrum.
% (version 1.0)
%
% "data" must be a 1-D signal.
% "dt" must be the time interval of signal
% "which_spectrum"
%       if you put 'a' the output is just amplitude spectrum
%       if you put 'p' the output is just phase spectrum
%       if you put 'ap' the output is both amplitude and phase spectrums (default)
% "color" you can specify the color of plots for example 'black', 'red', ...
%     default is 'blue' (good option for comparison of two signals)
% "threshold" you can choose a threshold value (it must be positive) for 
%      removing some frequencies with amplitudes lower than the threshold value; 
%      its default value is zero
% "deg_OR_rad"
%           if you put 'deg' the angle of phase spectrum is based on degree (default)
%           if you put 'rad' the angle of phase spectrum is based on radian
%
%-------------------------------Example------------------------------------
%
% n=0:1:999;
% dt=0.001;
% signal=cos(2*pi*10*n*dt+50)+cos(2*pi*73*n*dt-121)+cos(2*pi*184*n*dt+30)+...
%       +cos(2*pi*218*n*dt-57)+cos(2*pi*346*n*dt+74);
% 
% ft_spect(signal,0.001,'ap','red',0.1,'rad')
%

Answer 3

The following answer is based on Serge's solution; however, I want to add a bit more. The phase plot is very sensitive to the data; any minute change to the time vector or little introduction of noise will produce an entirely different phase plot. The phase plot in the question seems to show that the time vector is not exactly perfect. I include the following code to illustrate my point.

Fs=100; %sampling frequency
T=1; %signal length
N=T*Fs;     %number of samples
t=linspace(0,T-1/Fs,N);% do not use t=linspace(0,T,N) eventhough it is almost the same
f=-Fs/2:Fs/N:Fs/2-Fs/N; %frequency vector
x=cos(2*pi*10*t);
Xf=fft(x);
Xf=fftshift(Xf);
subplot(3,1,1)
plot(t,x) %original signal
subplot(3,1,2)
plot(f,abs(Xf)) %magnitude
subplot(3,1,3)
Xf(abs(Xf)<1e-6)=0; % remove tiny noise in Xf, crucial for plotting phase
plot(f,angle(Xf)) %phase

Notice in the phase plot that there are two small peaks that is only about 2e-15. They are just errors and should be zeros because the phase plot reflects the phase shift of a cosine function. There is no phase shift in this case, so the phase should be zeros across all frequencies.