fourier transform

[PG1995]

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.

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

 

[Ratchit]

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.

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

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

 

[steveB]

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.

 

[PG1995]

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.

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

??? 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

 

[steveB]

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.

https://mathworld.wolfram.com/SimpsonsRule.html
https://en.wikipedia.org/wiki/Simpson's_rule
https://en.wikipedia.org/wiki/Trapezoidal_rule

 

[PG1995]

Hi,

The Cosine Fourier Transform,

[latex]F(\omega )=\frac{2}{\pi }\int_{0}^{\infty }f(x)\cos (\omega x)d\omega[/latex]

of f(x)=1/√x is: √(2/πω)=√(/π²f). (Source: http://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.

 

[steveB]

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.

 

[MrAl]

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.

 

[PG1995]

Something seems wrong here. What is the variable "n" representing.

That's pi.

Also, the integral should be over the dx, not dw.

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

Cosine Fourier Transform:
[latex]F(\omega )=\frac{2}{\pi }\int_{0}^{\infty }f(x)\cos (\omega x)dx[/latex]
Inverse Cosine Fourier Transform:
[latex]f(x)=\int_{0}^{\infty }F(\omega )\cos (\omega x)d\omega[/latex]

So, again, the cosine fourier transform of f(x)=1/√x is: √(2/πω)=√(/π²f). (Source: http://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.

 

[MrAl]

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:
[latex]F(\omega )=\frac{2}{\pi }\int_{0}^{\infty }f(x)\cos (\omega x)dx[/latex]
Inverse Cosine Fourier Transform:
[latex]f(x)=\int_{0}^{\infty }F(\omega )\cos (\omega x)d\omega[/latex]

So, again, the cosine fourier transform of f(x)=1/√x is: √(2/πω)=√(/π²f). (Source: http://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.

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.

 

[steveB]

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.

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.

 

[steveB]

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.

 

[PG1995]

Thank you, MrAl.

The second expression makes no sense.

Are you talking about inverse fourier transform expression?

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.

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

 

[MrAl]

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

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.

 

[PG1995]

Thank you, Steve.

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.

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

 

[steveB]

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.

 

[PG1995]

Thank you, Steve.

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.

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

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

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

 

[MrAl]

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.

 

[steveB]

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").

 

[PG1995]

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.