How to plot phase and amplitude spectrum after doing fourier transform?

Hello, I am a new MATLAB user. I had a function which I did Fourier Transform for, and the result was: X(w)=1/(1+jw) where w is the frequency and " j " is the known imaginary number. I would like to know what code I should input in MATLAB in order to plot the phase and amplitude spectra of X(w). Thanks in advance...

Answer 1

[sixwwwwww]

Dear Bilal, You can plot in the following way the amplitude and phase of X:

figure, plot(w, abs(X)), title('Amplitude plot')
figure, plot(w, angle(X)), title('Phase plot')

Good luck!

[Bilal]

Hello and Thanks, however I am still getting an error! Here is my code:

>> w=sym('w'); >> X=(1/1+i*w); >> figure, plot(w, abs(X)), label('Amplitude plot') figure, plot(w, angle(X)), label('Phase plot')

The error is saying cannot convert to double from sym! Please help I am new to this

[sixwwwwww]

Do you have values for frequency? because you need some frequency values to plot amplitude and phase. In the following code I am using frequency values from 0-100:

syms w 
X = 1 / (1 + 1j * w); 
w_values = 0:100;
X_values = double(subs(X, w, w_values));
figure, plot(w_values, abs(X_values)), title('Amplitude plot'), xlabel('Frequency'), ylabel('Amplitude') 
figure, plot(w_values, angle(X_values)), title('Phase plot'), xlabel('Frequency'), ylabel('Phase')

I hope it is what you needed. Further more you can vary frequency values as you desire. Good luck!

[Bilal]

THANK YOU ! Well appreciated

[sixwwwwww]

You are welcome. Good luck!

[Sergio E. Obando]

Unless X is coming from a more complicated symbolic manipulation, you can plot with just:

w = 0:100;
X = 1 ./ (1 + 1j * w);
figure, plot(w, abs(X)), title('Amplitude plot'), xlabel('Frequency'), ylabel('Amplitude') 

figure, plot(w, angle(X)), title('Phase plot'), xlabel('Frequency'), ylabel('Phase')

Answer 2

[lakshya wadhwa]

syms w

X = 1 / (1 + 1j * w);

w_values = 0:100;

X_values = double(subs(X, w, w_values));

figure, plot(w_values, abs(X_values)), title('Amplitude plot'), xlabel('Frequency'), ylabel('Amplitude')

figure, plot(w_values, angle(X_values)), title('Phase plot'), xlabel('Frequency'), ylabel('Phase')

Answer 3

[Afshin Aghayan]

look at the following Matlab function, it can calculate phase spectrum as well as amplitude spectrum with a perfect accuracy:

https://www.mathworks.com/matlabcentral/fileexchange/63965-amplitude-and-phase-spectra-of-a-signal--fourier-transform-

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 most signals that I have tried, 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, the code filters result of FT by thresholding, it is very important step! 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 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

Answer 4

[rahul jain]

hey there,

i am new to matlab and i am stuck at this code that i have to do for research in optical systems.

i have quite a bit of code already which lets me plotthe zernike polynomial and its poin spread function(individually).But what im trying to do is to find a way to plot all of the zernike polynials in a single figure, and in the next figure their PSFs. also how i do i plot the phase of a zernike wavefont.

i will attach my code below, and hopefully someone (perhaps you) knows it better that i do.

any help and recommendations would be deeply appriciated.

close all; clear variables; clc;

pixels=512*1;

d=5000; %eye pupil size in um
feye=22200; %eye axial length (schematic eye model)
neye=1.33; %eye refractive index
lambda=0.55; %wavelength of light in um

%ZERNIKE COEFFICIENTS (organised according to the OSA/ANSI standard)
c1=0; %tilt(in um)
c2=0; %tip (in um)
c3=0; %oblique astigmatism (in um)
c4=0.5; %defocus (in um)
c5=0; %vertical astigmatism (in um) 
c6=0; %vertical trefoil (in um)
c7=0; %vertical coma (in um)
c8=0; %horizontal coma (in um)
c9=0; %oblique trefoil (in um)
c10=0; %oblique quadrafoil (in um)
c11=0; %oblique secondary astigmatism (in um)
c12=0; %primary spherical (in um)
c13=0; %vertical secondary astigmatism (in um)
c14=0; %vertical quadrafoil (in um)
cvector=[c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14]

%generated rms of wavefront
rms=0;
for ii=1:14
    rms=(cvector(ii))^2 + rms;
end
rms=sqrt(rms)

L1=10000; %pupil image size in um
L2=lambda*pixels*(feye/neye)/L1  %retinal image size in um 

wavefront=zeros(pixels,pixels); %Zernike wavefront
pupil=zeros(pixels,pixels); %field in pupil plane with aberrations
refpupil=zeros(pixels,pixels); %reffield in pupil plane without aberrations
psf=zeros(pixels,pixels); %Intensity PSF
reconstructwavefront=zeros(pixels,pixels); %least-square wavefront reconstruction

%calculating the field in the pupil plane
for ii=1:pixels
    yy=L1*(ii-pixels/2 -0.5)/pixels;
    for jj=1:pixels
        xx=L1*(jj-pixels/2 -0.5)/pixels;
        rho=sqrt(xx^2 + yy^2);
        if rho<=(d/2)
            rho0=rho/(d/2);
            theta=atan2(yy,xx);
                
            wavefront(ii,jj)=c1*2*rho0*sin(theta);
            wavefront(ii,jj)=c2*2*rho0*cos(theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c3*sqrt(6)*((rho0^2)*sin(2*theta))+wavefront(ii,jj);
            wavefront(ii,jj)=c4*sqrt(3)*(2*(rho0^2) -1)+wavefront(ii,jj);
            wavefront(ii,jj)=c5*sqrt(6)*((rho0^2)*cos(2*theta))+wavefront(ii,jj);
            wavefront(ii,jj)=c6*sqrt(8)*(rho0^3)*sin(3*theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c7*sqrt(8)*(3*(rho0^3)-2*rho0)*sin(theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c8*sqrt(8)*(3*(rho0^3)-2*rho0)*cos(theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c9*sqrt(8)*(rho0^3)*cos(3*theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c10*sqrt(10)*(rho0^4)*sin(4*theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c11*sqrt(10)*(4*(rho0^4)-3*(rho0^2))*sin(2*theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c12*sqrt(5)*(6*(rho0^4)-6*(rho0^2)+1)+wavefront(ii,jj);
            wavefront(ii,jj)=c13*sqrt(10)*(4*(rho0^4)-3*(rho0^2))*cos(2*theta)+wavefront(ii,jj);
            wavefront(ii,jj)=c14*sqrt(10)*(rho0^4)*cos(4*theta)+wavefront(ii,jj);
            
            pupil(ii,jj)=1*exp(1i*(2*pi/lambda)*wavefront(ii,jj)); %converting wavefront to phase
            refpupil(ii,jj)=1; %unaberrated reference pupil
            
        end
    end
end
dcut=round(pixels*d/L1); %size of non-zero pupil
wf=wavefront(round((pixels-dcut)/2 +1):round((pixels+dcut)/2),round((pixels-dcut)/2 +1):round((pixels+dcut)/2)); %image cut to match pupil
figure(1), imagesc(flipud(wf)), axis square, title('Wavefront [um]'), axis off, colormap('jet'), colorbar;

sensor=rot90(fftshift(fft2(pupil)),2);  %field in image sensor plane viewed head-on

disp(sensor)
intensity=(abs(sensor).^2)/max(max(abs(fftshift(fft2(refpupil)).^2)));  %this is the ratio of the aberrated PSF / unaberrated PSF in the image plane
figure(2), imagesc(intensity), axis square, title(append('PSF viewed head-on',' ',mat2str(round(L2)),'x',mat2str(round(L2)),' um2')), axis off, colormap('gray'), colorbar;

Answer 5

[Siri]

Plot the spectra (s) as a function of wavelength (lambda). Use point markers (.) and a solid line (-) to connect the points.

Add the x-label "Wavelength" and the y-label "Intensity" to the plot.