1. Welcome to our site! Electro Tech is an online community (with over 170,000 members) who enjoy talking about and building electronic circuits, projects and gadgets. To participate you need to register. Registration is free. Click here to register now.
    Dismiss Notice

fourier transform

Discussion in 'Mathematics and Physics' started by PG1995, Nov 8, 2017.

  1. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    Hi,

    Could you please help me to fix the code below? I know that the code is completely wrong but I have tried to write it down so that you can get an idea what I'm trying to do.

    I'm trying to understand the fourier transform and I needed the code to formulate my questions. Thank you.

    Code (text):

        %%file name fourier.m
        clear all; close all; clc;
        f = 0; %% 0 Hz to 100 Hz, with 0.1 step, there would be 1000 frequencies
        cos_v=[1:1:1000]; %%vector to store values
        f_w=[]; %% defining a vector

        for i1=[0:1000]; %% for loops begins, with loop index i1
        f_w = sqrt(1/(4*f));
        cos_v(i1)=f_w*cos(2*pi*f*x); %% i1 serves as an index for cos_v vector elements too
        f=f+0.1;
        end
     
        combined_sinusoids=sum(cos_v); %% cos_v(1)+cos_v(2)+cos_v(3)+...
        x = [0:0.2:500];
        y = combined_sinusoids;
        plot(x, y)
     
     
    Last edited: Nov 8, 2017
  2. Ratchit

    Ratchit Well-Known Member

    Joined:
    Mar 12, 2008
    Messages:
    2,012
    Likes:
    93
    At least you could disclose what language the code is written. Do you expect us to get up to speed on that particular language and debug your code? Perhaps a more direct question would be appropriate.

    Ratch
     
  3. steveB

    steveB Well-Known Member Most Helpful Member

    Joined:
    Jan 16, 2009
    Messages:
    1,307
    Likes:
    639
    I see this is coded in Matlab, and there will be many people here that are familiar with Matlab coding, but I have to echo the basic sentiment that you haven't provided us a good effort in asking this question. When I run this code, I get an error on line 9 that says x is not defined. I see it is defined later in the code, but why not before? Basically, it is not clear at all what you are trying to do. It would be better for you to provide the equations you wrote out before you started coding. ... Hmm, did you do that? If not, that's a big mistake. Your method should be well planned out before you even begin coding.

    That said, what definition of "Fourier Transform" are you using? It looks like you are only considering the cosine component and not the sine component. The usual Fourier transform is going to involve complex values from the i*sin piece, or the exp (i x) if you do it from that definition.

    To reiterate, provide the definition of the Fourier transform (continuous or discrete, exp or sin/cos etc.), then a derivation of how you get to the form you will implement in code (i.e integrals converted into summations , variable definitions, index definitions etc.). And, most importantly, the code you do provide should run without errors.
     
    Last edited: Nov 9, 2017
    • Like Like x 1
  4. dave miyares

    Dave New Member

    Joined:
    Jan 12, 1997
    Messages:
    2
    Likes:
    -10


     
  5. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    Thank you, steveB.

    Ratchit: It's Matlab code. Sorry, I should have mentioned it.

    You could check this PDF to see what I'm trying to do. I have been trying to get a conceptual understanding of Fouier transform like how it basically works and why it works as a kind of series representation of a function.

    How could I transform FORM #2 into FORM #1? Thank you.

    Re: Code
    I believe MATLAB treats every variable (even a constant) as a vector. A vector is a one-dimensional array and a matrix is two-dimensional array. I've experience with C++ so I'd tend to stick twith he syntax of C++ so as I learn more about MATLAB I'll start thinking more in MATLAB terms.

    Given below is the code with some revisions; also the error.

    Code (Matlab):

    %%file name fourier.m
        clear all; close all; clc;
        f = [0.1]; %%defining and initializing a constant vector
                   %%0 Hz to 100 Hz, with 0.1 step, there would be 1000 frequencies
        cos_v=zeros(1000); %%defining and initializing vector to store values
                           %%a vector with 1000 zeroes
        f_w=[0]; %% defining and initializing a vector, fourier transfor,
        syms x; %%"x" is a floating point variable which can take values from 0 to
                %%infinity, in C++ we could have written "float x"

        for i1=[1:1000] %% for loops begins, with loop index i1
        f_w = sqrt(1/(4*f)); %% fourier transform of 1/sqrt(x)
        cos_v(i1)=f_w.*cos(2*pi*f*x); %% i1 serves as an index for cos_v vector elements too
        f=f+0.1; %%incrementing f
        f_w=0; %%setting f_w=0 for next iteration
        end
       
        combined_sinusoids=sum(cos_v); %% cos_v(1)+cos_v(2)+cos_v(3)+...
        y = combined_sinusoids;
       
        x = input('Enter value of x: ');
       
        fprintf('The value of combined_sinusoids is %d \n',y)
       
        fprintf('The value of 1/sqrt(x) is %d \n',1/sqrt(x))
       
        fprintf('Difference between 1/sqrt(x) and approximation by combined_sinusoids is %d \n',abs(y-(1/sqrt(s))))
     
    Code (Matlab):

    ??? The following error occurred converting from sym to double:
    Error using ==> mupadmex
    Error in MuPAD command: DOUBLE cannot convert the input expression into a double array.

    If the input expression contains a symbolic variable, use the VPA function instead.

    Error in ==> fourier_new at 13
        cos_v(i1)=f_w*cos(2*pi*f*x); %% i1 serves as an index for cos_v vector elements too

     
     
    Last edited: Nov 12, 2017
  6. steveB

    steveB Well-Known Member Most Helpful Member

    Joined:
    Jan 16, 2009
    Messages:
    1,307
    Likes:
    639
    It's still not 100 percent clear what you are trying to do, but I will take a guess and you can tell me if I'm correct. Before that, I will mention that your code is trying to use symbolic math commands. I don't have the symbolic toolbox installed, so I can't run the code, but I wonder why you are using symbolic processing. If the symbolic processor can do it, then you can do it too, by hand. Usually Fourier transform programs are intended to be numerical to handle general cases that can't be solved symbolically. Also, you are using a Cosine Fourier Transform and not the usual Fourier Transform. I'm not sure why you are using that one, but it is fine because it is just an integral that is define for functions that allow it to converge.

    OK, so I think what you are trying to do is simply calculate a numerical estimation of a one dimensional integral. This should be straightforward if you convert the integral to an approximate summation using any of the standard numerical integral techniques. Have you ever used the "Trapazoidal Rule" or "Simpson's Rule"? I would recommend this approach, and start with the trapazoidal rule, which is simpler. Then just code it in a loop. No need for anything complicated.

    http://mathworld.wolfram.com/SimpsonsRule.html
    https://en.wikipedia.org/wiki/Simpson's_rule
    https://en.wikipedia.org/wiki/Trapezoidal_rule
     
    • Like Like x 1
  7. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    Hi,

    The Cosine Fourier Transform,



    of f(x)=1/√x is: √(2/πω)=√(/π²f). (Source: https://www.math.uni-sb.de/ag/fuchs/PDE14-15/pde14-15-lecture-7.pdf)

    When f=0.1 Hz, Fourier transform is √(/π²f)=1.007. What does it really mean? In my opinion, it means that the magnitude of cosine sinusoid with frequency 0.1 Hz which combines with other cosine sinusoids to make up the signal f(x)=1/√x is 1.007. Do I have it correct? Thank you.
     
  8. dave miyares

    Dave New Member

    Joined:
    Jan 12, 1997
    Messages:
    2
    Likes:
    -10


     
  9. steveB

    steveB Well-Known Member Most Helpful Member

    Joined:
    Jan 16, 2009
    Messages:
    1,307
    Likes:
    639
    Something seems wrong here. What is the variable "n" representing. Also, the integral should be over the dx, not dw.

    It would be good to show the Fourier transform and its proper inverse integral as well. The inverse integral will answer your question. You basically have it correct, but there may be an additional constant involved, which would become apparent in the inverse transform.
     
    • Like Like x 1
  10. MrAl

    MrAl Well-Known Member Most Helpful Member

    Joined:
    Sep 7, 2008
    Messages:
    11,077
    Likes:
    964
    Location:
    NJ
    Hi,

    That's pi.

    There are a lot of uses for the cos transform but i cant remember many of them.

    One use is for encoding the jpeg file format. The idea is to transform the two dimensional picture data into frequency components, and those components very more smoothly than the raw picture data, so when the next stage of the jpeg decoding comes up which uses differential coding, the new data changes much less and so the numerical data is actually smaller and codes with less bits. That makes the whole file much smaller even with a quality factor of 100 percent which is lossless.

    The way to understand the whole process of using the cos transform though is to look at using both the transform and the inverse transform for just one or a few values of 'x' or time 't'. When you look at the reconstruction (the inverse transform) you see the components adding back up to equal the original time (or distance, or pixel) data.

    I have not checked your actual results yet though or that site you link to.
     
    Last edited: Nov 20, 2017
    • Like Like x 1
  11. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    That's pi.

    You're right. It should have been dx. The mentioned source document had it wrong and I didn't notice it.

    Cosine Fourier Transform:

    Inverse Cosine Fourier Transform:


    So, again, the cosine fourier transform of f(x)=1/√x is: √(2/πω)=√(/π²f). (Source: https://www.math.uni-sb.de/ag/fuchs/PDE14-15/pde14-15-lecture-7.pdf)

    When f=0.1 Hz, Fourier transform is √(/π²*f)=1.007. What does it really mean? In my opinion, it means that the magnitude of cosine sinusoid with frequency 0.1 Hz which combines with other cosine sinusoids to make up the signal f(x)=1/√x is 1.007. Do I have it correct? Thank you.
     
  12. MrAl

    MrAl Well-Known Member Most Helpful Member

    Joined:
    Sep 7, 2008
    Messages:
    11,077
    Likes:
    964
    Location:
    NJ
    The second expression makes no sense.
    Also, when you take the square root of pi squared you get pi:
    1/(pi*sqrt(f))


    the expression /pi^2*f does not make sense on several levels.
     
    • Like Like x 1
  13. steveB

    steveB Well-Known Member Most Helpful Member

    Joined:
    Jan 16, 2009
    Messages:
    1,307
    Likes:
    639
    So, I think you have the basic concept correct, but I hesitate to say your are correct. The inverse transform is the key to your answer, and that's why I wanted you to state it. The inverse transform is saying that the integration of the Fourier-cosine transform reconstructs the signal. If this were a summation, rather than an integral, I would say your interpretation is correct. With the integral, you need to be more careful in a proper statement. But, if you are just trying to get a mental picture of what is happening, then you can rewrite the inverse integral transform as an approximate summation transform. The dw gets replaced by a delta-w and the integration limits get replaced by an integer summation index (e.g. integer k=0, 1, 2, ... infinity) and the w gets replaced k*delta-w. In this form, your interpretation can be that the k'th frequency component (for the cosine) has a magnitude F(k*delta-w)*delta-w.
     
    Last edited: Nov 21, 2017
    • Like Like x 1
  14. steveB

    steveB Well-Known Member Most Helpful Member

    Joined:
    Jan 16, 2009
    Messages:
    1,307
    Likes:
    639
    By the way, you may want to look up "completeness" for the Fourier-cosine transform. I have a vague recollection that this is not a complete transform and may only be complete for "even" functions. Normally, the full Fourier transform with has both sine and cosine components is considered to be complete.
     
    • Like Like x 1
  15. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    Thank you, MrAl.

    Are you talking about inverse fourier transform expression?

    Sorry, I meant to write cosine fourier transform of f(x)=1/√x is: √(2/πω)=√(1/π²f)=1/π*√(1/f).
     
  16. MrAl

    MrAl Well-Known Member Most Helpful Member

    Joined:
    Sep 7, 2008
    Messages:
    11,077
    Likes:
    964
    Location:
    NJ
    Hi,

    No need to be sorry, but if you want people to know exactly what you want to ask then you at the very very least have to check over your syntax.

    The inverse cosine transform. See how it combines back together.

    This might be easier to see though if you look at the discrete cosine transform.
     
    • Thanks Thanks x 1
  17. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    Thank you, Steve.

    I think that you are right. In an online resource it was written that fourier transform reduces to just cosine part, i.e. cosine fourier transform, for even functions. So, I was also wondering that what if the function is not even. For example, the function f(x)=1/sqrt(x) being discussed is not an even function so this was already bothering me that how come it's cosine fourier transform is an exact representation of 1/sqrt(x) in frequency domain (btw 1/sqrt(|x|) is an even function).
     
  18. steveB

    steveB Well-Known Member Most Helpful Member

    Joined:
    Jan 16, 2009
    Messages:
    1,307
    Likes:
    639
    Perhaps one of the things that protects you here is that the integration limits are only on the positive part of the x-axis and not on the negative part. This converts your function (effectively) into an even function by ignoring the part in the negative domain. If your function is not even, or if the function is undefined in the negative domain, then you never see the problem from this fact.

    Anyway, I'm not a mathematician and I only wanted to point this out so that you can be wary in the future.
     
    • Like Like x 1
  19. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    Thank you, Steve.

    But I don't think that you could reconstruct the function f(x) using its fourier transform. MrAl already noted that the inverse fourier transform above doesn't make any sense. I tried to write a Matlab code to test it using 50 million sinusoids and it looks like that f(x)=1/sqrt(x) cannot be reconstructed from its cosine fourier transform F(f)=(1/pi)*sqrt(1/f).

    Code (Matlab):

    %file name fourier.m

          clear all; close all; clc;
          f = [0.001]; %defining and initializing a constant vector
                   %0 Hz to 100 Hz, with 0.1 step, there would be 1000 frequencies
          cos_v=0.*[1:50000000]'; %defining and initializing vector to store values
                           %a vector with 1000 zeroes
          f_w=[0]; % defining and initializing a vector, fourier transform,
             
          user_decision = 'y';
       
    while user_decision == 'y'
        x = input('Enter value of x: ');
     
        for i1=[1:50000000] % for loops begins, with loop index i1
          f_w = (1/pi)*sqrt(1/f); % fourier transform of 1/sqrt(x)
          cos_v(i1)=f_w.*(cos(2*pi*f*x)); % i1 serves as an index for cos_v vector elements too
          f=f+0.001; %incrementing f
          f_w=0; %setting f_w=0 for next iteration
        end
     
        combined_sinusoids=sum(cos_v); % cos_v(1)+cos_v(2)+cos_v(3)+...
        y = combined_sinusoids;
     
        fprintf('The value of combined_sinusoids is %d \n',y)
     
        fprintf('The value of 1/sqrt(x) is %d \n',(1/sqrt(x)))
     
        fprintf('Difference between 1/sqrt(x) and approximation by combined_sinusoids is %d \n',abs(y-(1/sqrt(x))))
     
        user_decision=input('Would you like to continue, enter Y/N: ','s');
     
        if user_decision=='n'
            break
        end

    end
     
    For x=2, the result is:
    Code (text):

    Enter value of x: 2
    The value of combined_sinusoids is 9.784050e+01
    The value of 1/sqrt(x) is 7.071068e-01
    Difference between 1/sqrt(x) and approximation by combined_sinusoids is 9.713339e+01
    Would you like to continue, enter Y/N: n
     
     
  20. MrAl

    MrAl Well-Known Member Most Helpful Member

    Joined:
    Sep 7, 2008
    Messages:
    11,077
    Likes:
    964
    Location:
    NJ
    Hi,

    When i said something did not make sense it was the "/denom" that did not make any sense but you cleared that up stating that it was really supposed to be "1/denom" which does make sense. Whether it is correct or not i did not look into.
     
    • Thanks Thanks x 1
  21. steveB

    steveB Well-Known Member Most Helpful Member

    Joined:
    Jan 16, 2009
    Messages:
    1,307
    Likes:
    639
    PG,

    The way you can check whether the function can be reconstructed is to use the completeness theorem for the Proper Fourier Transform. For example, you can take your 1/sqrt(x) function and turn it into an even function by making a mirror image about the y-axis. Then you can apply the Fourier transform to that. If you find that the Fourier transform has only a value from the "cosine part" and a zero from the "sine part", then I think you have a reconstruction possible with the Cosine Fourier inverse integral. Note that you will reconstruct the new (mirror imaged) function, and not the original function. In other words, if the Sine Fourier transform does not universally give you zero, the Cosine transform will not allow you to reconstruct the signal because it needs the Sine part in order to be complete.

    Keep in mind that doing a numerical reconstruction is a procedure requiring great care and error bounding. Otherwise, if you find you can't reconstruction the signal, you never know if it is due to a coding mistake, or accumulation of numerical errors, or a theoretical problem with the method (i.e. the classic logical problem with trying to prove a "negative").
     
    • Like Like x 1
  22. PG1995

    PG1995 Active Member

    Joined:
    Apr 18, 2011
    Messages:
    1,681
    Likes:
    13
    Hi,

    Taylor series is a kind of power series where coefficients are derivatives. A function f(x) is approximated using Taylor series around a fixed point, 'a', and as the function is evaluated at points away from 'a', error would increase. I'm assuming finite number of terms. In other words, Taylor series approximates a function locally around a fixed point. On the other hand, Fourier series or transform approximates a function over its defined domain and not around a fixed point. In other words, Fourier series or transform approximation is more of global. I understand comparing Taylor series and Fourier series is more like comparing oranges and apples but I'm trying to make a very general comparison between the two.

    Taylor series was introduced around 1715 (https://en.wikipedia.org/wiki/Taylor_series). Calculus was developed around 1680. The following is an excerpt from https://en.wikipedia.org/wiki/Calculus#Modern, "In other work, he developed series expansions for functions, including fractional and irrational powers, and it was clear that he understood the principles of the Taylor series. He did not publish all these discoveries, and at this time infinitesimal methods were still considered disreputable." Does the part in red mean that Newton had experimented with a form of Taylor series in his work even before the series was formally introduced by Brook Taylor?

    Thank you for your help.
     

Share This Page