Approximating the inverse FT of a unit pulse using a Riemann sum

[PainterGuy]

TL;DR Summary

I have been trying to write a MATLAB code to approximate an inverse trigonometric Fourier transform of unit pulse using Riemann sum. My code isn't working.

Hi,

Although I'm using trigonometric form of Fourier transform, first I'd discuss both, exponential and trigonometric forms, for the sake of context.

Now proceeding toward the main question and we would only be using trigonometric form.

% file name unit_pulse_fourier_transform_riemann_sum.m
% riemann sum: https://www.math.hmc.edu/calculus/tutorials/riemann_sums/gif/figure3.gif

clear all; close all; clc;

Dw=0.5; %delta Omega
i=1; %frequency increases 0 to 500 in steps of 0.5 Hz therefore max i should be 1000;
w=0.*[0:0.5:500];
w_mid=zeros(1000);
x=linspace(-2*pi,2*pi,13000); % dividing the interval into 13000 points
Rs=zeros(1000);
Rs_x=zeros(13000);

for i1=1:13000

    w(1)=0;
    w(2)=0.5;

    for i=1:1000
   
    w_mid(i)=(w(i+1) + w(i))/2 %w_mid(1) is first middle point

    Rs(i)=(2*sin(pi*w_mid(i))/w_mid(i))*cos(x(i1)*w_mid(i))*Dw;

    w(i+2)=w(i+1) + Dw;
   
    end

Rs_x(i1)=sum(Rs); % Rs_x holds values for evaluated Rs

end

plot(x,Rs_x);
hold;

I understand that the code above involves many iterations so to simplify it just to see if it works, I reduced the number of iterations. It doesn't work. I have include the errors at the bottom. Could you please help me to fix the code? Thank you!

% file name unit_pulse_fourier_transform_riemann_sum.m
% for riemann sum: https://www.math.hmc.edu/calculus/tutorials/riemann_sums/gif/figure3.gif

clear all; close all; clc;

Dw=0.5;
w=0.*[0:0.5:25];
w_mid=zeros(50);
x=linspace(-2*pi,2*pi,60);
Rs=zeros(50);
Rs_x=zeros(60);

for i1=1:60

    w(1)=0;
    w(2)=0.5;

    for i=1:50
   
    w_mid(i)=(w(i+1) + w(i))/2 %w_mid(1) is first middle point

    Rs(i)=(2*sin(pi*w_mid(i))/w_mid(i))*cos(x(i1)*w_mid(i))*Dw;

    w(i+2)=w(i+1) + Dw;
   
    end
   
Rs_x(i1)=sum(Rs);

end

plot(x,Rs_x);
hold;

In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in delete_fourier_pulse (line 28)
Rs_x(i1)=sum(Rs);

 

[RPinPA]

Do you have any experience in debugging Matlab?

It's hard to diagnose this from looking at the code. Add the command "dbstop if error" after line 4, since your "clear all" clears all breakpoints and they need to be reset.

If you get the error, check the different parts of that assignment, such as the value of i1, the value of sum(Rs), the size of Rs_x, etc. See if they are all what you expect.

Another option is to embed the problem assignment statement in a TRY...CATCH...END block and then put some diagnostics in the CATCH block.

 

[PainterGuy]

Thank you for your attempt to help me!

I didn't apply any of the debugging and don't have any experience with it either. Sadly, I only use MATLAB when I don't see any other way and only use it quite rarely.

I was able to trace the problem. The main problem was with the lines such as "w_mid=zeros(50);". I was under the impression that it creates a row vector with zero values but it actually generates a matrix.

Also, we also need to use the factor of 1/π and in my first post I did't even mention it. I have updated the code. Without the factor of 1/π, the pulse height was close to "3" when it should have been just "1".

% file name unit_pulse_fourier_transform_riemann_sum_reduced_dimensions.m
% for riemann sum: https://www.math.hmc.edu/calculus/tutorials/riemann_sums/gif/figure3.gif

clear all; close all; clc;

Dw=0.5;
w=0.*[0:0.5:25];
x=linspace(-2*pi,2*pi,60);
Rs_x=0.*[1:60];

for i1=1:60

    w_mid=0.*[1:50];
    Rs=0.*[1:50];
    w(1)=0;
    w(2)=0.5;

    for i=1:50

    w_mid(i)=(w(i+1) + w(i))/2; %w_mid(1) is first middle point

    Rs(i)=(1/pi)*(2*sin(pi*w_mid(i))/w_mid(i))*cos(x(i1)*w_mid(i))*Dw;

    w(i+2)=w(i+1) + Dw;

    end

Rs_x(i1)=sum(Rs);

end

plot(x,Rs_x);
hold;

I have a question about MATLAB in general. I believe all of the calculations in MATLAB are based on discrete values. I mean that you need to input discrete values because MATLAB can only work with matrices. In other words, you always need to discretize any problem you are trying to solve using MATLAB. Do I make sense?

Thanks a lot!Note to self:
If you try to plot Fourier transform of unit pulse, i.e. 2sin(ωπ)/ω, you need to create an exception to handle ω=0 case because unless the limit is taken, the expression cannot be evaluated. See the example code below. There is no goto statement in MATLAB.

% file name use_of_continue_statement.m

clear all; close all; clc;

w=linspace(-2*pi,2*pi,40);

for i=1:40
   
    if w(i)==0
        FT(i)=1;
        continue
    end
  
    FT(i)=(2*sin(w(i)*pi))/w(i);

end

plot(w,FT);
hold;

Helpful links:
1: http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/ref/continue.html
2: https://www.mathworks.com/help/matlab/ref/break.html
3: https://stackoverflow.com/questions/1082605/jump-command-in-matlab

 

[RPinPA]

I was able to trace the problem. The main problem was with the lines such as "w_mid=zeros(50);". I was under the impression that it creates a row vector with zero values but it actually generates a matrix.

I should have caught that as I have extensive experience with Matlab and I'm well aware that zeros(50) means a 50x50 matrix. But my normal approach to debugging is to set a breakpoint, then try to build the problem expression bit by bit and see where it fails, and then use that as a clue to what is causing the failure.

Also, we also need to use the factor of 1/π and in my first post I did't even mention it. I have updated the code. Without the factor of 1/π, the pulse height was close to "3" when it should have been just "1".

Again I was well aware there were normalization constants left out but like you considered it unimportant since it wasn't the source of the error. There are many different forms of the FT with different normalizations and I rarely worry about the scaling when doing a transform.

I have a question about MATLAB in general. I believe all of the calculations in MATLAB are based on discrete values. I mean that you need to input discrete values because MATLAB can only work with matrices. In other words, you always need to discretize any problem you are trying to solve using MATLAB. Do I make sense?

Yes you make perfect sense. The answer is that MATLAB does include a symbolic toolbox as an option, which I think is built on MAPLE. It can work on algebraic expressions, do symbolic differentiation and integration, etc.
https://www.mathworks.com/products/symbolic.html
Toolboxes are add-ons, so you have to check if your installation includes this particular toolbox. Issuing the VER command will tell you the versions of Matlab and all installed toolboxes.

 

[PainterGuy]

Hi again,

I was trying to apply the approach of approximation using Riemaan sum to a shifted unit pulse by π/2. It looks like I'm doing the math right but there is a problem with the plot at bottom.

In the following code, I'm using frequencies from 0 Hz to 100 Hz using increments of 0.5 Hz, dividing the x-axis into 120 points from -4π to 4π. The rest of code is same as in my previous post.

% file name shifted_fourier_transform_unit_pulse_2.m
% for riemann sum: https://www.math.hmc.edu/calculus/tutorials/riemann_sums/gif/figure3.gif
% using shift of x=t=pi/2

clear all; close all; clc;

Dw=0.5;
w=0.*[0:0.5:100];
x=linspace(-4*pi,4*pi,120);
Rs_x=0.*[1:120];

for i1=1:120

    w_mid=0.*[1:200];
    Rs=0.*[1:200];
    w(1)=0;
    w(2)=0.5;

    for i=1:200

    w_mid(i)=(w(i+1) + w(i))/2; %w_mid(1) is first middle point

    Rs(i)=(1/pi)*((cos((pi/2)*w_mid(i))*(2*sin(pi*w_mid(i))/w_mid(i))*cos(x(i1)*w_mid(i)))+(sin((pi/2)*w_mid(i))*(2*sin(pi*w_mid(i))/w_mid(i))*(sin(x(i1)*w_mid(i)))))*Dw;

    w(i+2)=w(i+1) + Dw;

    end

Rs_x(i1)=sum(Rs);

end

plot(x,Rs_x);
hold;

You can see in the plot below that the unit pulse is okay but I don't understand why it falls down to "-1" in the yellow regions. It should have stayed zero everywhere except from -1.57 to 4.7 (or, -π + π/2 to π + π/2) on the x-axis. Could you please help me with that? Where am I having it wrong? Thank you!

 

[PainterGuy]

Hi again,

I was trying to fix the problem from my previous posting about the value falling to "-1". While doing this, I tried to see if the same problem happens with the code posted in post #3 for the un-shifted pulse. The same problem happens when I changed the limits for x. The plot is shown below.

Line #8 from the original code.

x=linspace(-2*pi,2*pi,60);

changed to

x=linspace(-4*pi,4*pi,60);

Then, I changed the following two lines from the same code of post #3 to see a larger picture of what was happening over an extended range of x values.

Line numbers #8 and #9 from the original code.

x=linspace(-2*pi,2*pi,60);
Rs_x=0.*[1:60];

changed to

x=linspace(-10*pi,10*pi,300);
Rs_x=0.*[1:300];

The plot below shows what I got. It gave rise to more questions. The plot seems to represent more of an inverse Fourier series! In Fourier transform, the assumption is made that the period is at infinity therefore there should have only been a single pulse and no other repetitions until we get to the infinity. In other words, it should have only been "1" from -π to π and everywhere else it should have been "0".

I'd really appreciate you could guide me where I'm going wrong. It's more about math now.

I also tried to see if something was wrong with the inverse Fourier transform integral,

, using numerical integration. It looks like integral is good.

By the way, in the code below, I'm using "quad" function to do the numerical integration and I understand that it's a deprecated function and "integral" function should be used instead. The "integral" function doesn't work for me and we can discuss it later.

I do get the following warning while using quad function but it worked nonetheless. The plot is shown below.

Warning: Maximum function count exceeded; singularity likely.
> In quad (line 98)
  In fourier_transform_unit_pulse_numerical_integration (line 18)

% file name "fourier_transform_unit_pulse_numerical_integration.m"
% using quad function

clear all; close all; clc;

x=linspace(-8*pi,8*pi,300);
Rs_x=0.*[1:300];for i1=1:300

    x1=x(i1);

    f=@(w)(1/pi)*(((2*sin(pi*w)/w)*cos(x1*w)));

    Rs_x(i1)=quad(f,0,60);

end

plot(x,Rs_x);
hold;

PS:
I believe that it has to do something with the value of Δω being used in the code for Riemann sum approximation; in the code Δω is represented as Dw. As the value of Dw is decreased, the period becomes larger; which would mean that as Dw goes to zero, the period becomes infinite. But I have not been able to understand how the value of Δω or Dw is affecting this behavior.

code: https://imagizer.imageshack.com/img921/8572/OaEQzV.jpg
plot: https://imagizer.imageshack.com/img923/5892/KYny42.jpg

code: https://imagizer.imageshack.com/img921/6286/2I99PH.jpg
plot: https://imagizer.imageshack.com/img924/4611/I99GMQ.jpg

 

[PainterGuy]

Hi again,

I have been trying to get over this confusion that why Riemann sum approximation to inverse trigonometric Fourier transform becomes inaccurate as the range of 'x' values is extended.

Toward the end of previous post, I had concluded that it had to do something with the value of Δω or Dw being used.

In the code below, "%%" means that the line has a variable which could be changed, or when you change a certain variable, try to go through these lines to see if any changes are needed. The three plots shown under the code were generated by changing value of Dw. In first plot, Dw=0.5, the first jump to "-1" occur around x=-10 and x=10. In second plot, Dw=0.25, the first jump to "-1" occur around x=-22 and x=22. In third plot, Dw=0.125, the first jump to "-1" occur around x=-47 and x=47. It clearly shows that it's the value of Dw which is causing this "jumping to -1" behavior. As the value of Dw becomes smaller, the jump occurs farther away.

% file name "unit_pulse_fourier_transform_riemann_sum_plotting_sinusoids.m"
% for riemann sum: https://www.math.hmc.edu/calculus/tutorials/riemann_sums/gif/figure3.gif
%{
In the code below inverse Fourier transform of unit pulse is found using
trigometric form of Fourier transform. The integral is approximated using Riemann sum.
Each sinusoid is plotted and then all the sinusoids are added to to plot an appeoximate pulse
%}

clear all; close all; clc;

Dw=0.5; %% you can change this value
max_freq_used=25; %% you can change this value
w=0.*[0:0.5:max_freq_used];
w_mid_max=max_freq_used/Dw;
w_mid=0.*[1:w_mid_max];
w(1)=0;
w(2)=Dw;
x=linspace(-20*pi,20*pi,500); %% you can change this value
y_values=zeros(w_mid_max,500); %%(max value of w_mid,number of total x intervals)for i1=1:w_mid_max   % max value of index is based on max w_midRs=0.*[1:500]; %%[1:number of x intervals]

w_mid(i1)=(w(i1+1) + w(i1))/2; %w_mid(1) is first middle point    for i=1:500 %% total number of x intervals

    Rs(i)=(1/pi)*(2*sin(pi*w_mid(i1))/w_mid(i1))*cos(x(i)*w_mid(i1))*Dw;
    y_values(i1,i)=Rs(i);

    end

plot(x,Rs);
hold on;
w(i1+2)=w(i1+1) + Dw;

end

z=sum(y_values);
plot(x,z,'linewidth',3);
hold;

Plot #1

Plot #2

Plot #3

Let's see why Dw is affecting this behavior. In Plot #3, the blue sinusoid is the one we start with and w_mid=(0+0.125)/2=0.0625. In other words, 2πf=0.0625 ⇒ f=0.00995 Hz ⇒ T=100.5. If we look closely at the blue sinusoid, it's period is actually 100 units. Think what would the period be if the frequency was, say, f=0.0000000000000000001 Hz?! The period would become almost infinite. In actual inverse Fourier transform, you have frequencies smaller than that! Also, if the period is, say, 100 units, it means for 50 units the sinusoid is going to be positive and for the next 50 units sinusoid is going to be negative. Therefore, when the period is infinite, these positive and negative half cycles would stretch over infinite number of units.

While using Riemann sum to approximate inverse Fourier transform, we increase the frequency using increments, for example in case of Dw=0.125, w_mid1=0.0625 (f=0.00995 Hz) then w_mid2=0.1875 (f=0.02984 Hz), then w_mid3=0.3125 (f=0.04974 Hz) and so on.

Let's look at the integral again

. In yellow we have a sinc function which gives the amplitude for cosine sinusoid whose frequency could vary from 0 to ∞. Remember in trigonometric form Fourier transform we only use positive frequencies. The code below could be used to plot sinc function. The plot is also shown. From the plot, we can see that the amplitude is maximum for ω cyclic frequencies between 0 < ω < 1. But in case of Riemann sum, we are only using increments which make us miss many 'important frequencies' between 0 < ω < 1. Further, think of a cos(ωx)=cos(2πf*x) where f=0.0000000000000000001 Hz,for such a sinusoid cosine is going to be positive or we can say its positive cycle is going to stretch over an infinite number of units. Using increments have us miss many such cosines whose positive half cycles extend over very, very large duration like a'positive' ceiling.

Note that for cos{2π(1)*x}, its positive half cycle extends over 0.5 units and its total period is 1 unit.

We can say that Riemann sum is good for integration and as a general approximate for area but when it is used to approximate inverse Fourier transform, it results into missing many important frequencies whose weightage or importance is way more compared to the other frequencies.

clc; clear all; close all;

w=[0:0.1:20];

for i=1:201

    if w(i)==0
    sinc(i)=2;
    continue
    end

sinc(i)=2*sin(w(i)*pi)/(w(i)*pi);

end

plot(w,sinc);
grid on;
hold;

Helpful link(s):
1: (www.youtube.com/watch?v=)9sd4DWragBg ; remove parentheses