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.
% 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
% 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
% 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