Computing Fast Fourier Transform in R












0












$begingroup$


This is a math question with programming application. So I am trying to find 3 things given a certain function in the x domain when transformed into the spectral domain.
1) the Amplitude
2) The Frequency
3) The Phase
In R (statistical software) I have coded the following function:



y=7*cos(2*pi*(seq(-50,50,by=.01)*(1/9))+32)
fty=fft(y,inverse=F)
angle=atan2(Im(fty), Re(fty))
x=which(abs(fty)[1:(length(fty)/2)]==max(abs(fty)[1:(length(fty)/2)]))
par(mfcol=c(2,1))
plot(seq(-50,50,by=.01),y,type="l",ylab = "Cosine Function")
plot(abs(fty),xlim=c(x-30,x+30),type="l",ylab="Spectral Density in hz")


I know I can compute the frequency by taking the bin value and diving it by the size of the interval(total time of the domain). Since the bins started at 1, when it should be zero, it would thus be
frequency=(BinValue-1)/MaxTime, which does get me the 1/9'th I have in the function above.



I have two quick questions:



First) I am having trouble computing the phase. the density function peaks at 12 (see bottom graph), shouldn't the value then the value of the phase be 2*pi+angle[12]
but I am getting a value of




angle[12]
[1] -2.558724




which puts the phase at 2*pi+angle[12]=3.724462. But that's wrong the phase should be 32 radians. What am I doing wrong?



Second) How do I convert the amplitude which is abs(fty) into the constant in my consine function above, which is 7? abs(fty)[12]=34351.41 , how do I convert this number to 7?. I appreciate your help!










share|cite|improve this question











$endgroup$












  • $begingroup$
    So, by $R$, you mean the statistical analysis package, and not the real numbers?
    $endgroup$
    – Thomas Andrews
    Sep 1 '14 at 17:37










  • $begingroup$
    Yes, the stat program. sorry about the confusion, i'll make the edit. thought this forum would be more apt than stack overflow because it really is more of a math question.
    $endgroup$
    – jessica
    Sep 1 '14 at 17:40
















0












$begingroup$


This is a math question with programming application. So I am trying to find 3 things given a certain function in the x domain when transformed into the spectral domain.
1) the Amplitude
2) The Frequency
3) The Phase
In R (statistical software) I have coded the following function:



y=7*cos(2*pi*(seq(-50,50,by=.01)*(1/9))+32)
fty=fft(y,inverse=F)
angle=atan2(Im(fty), Re(fty))
x=which(abs(fty)[1:(length(fty)/2)]==max(abs(fty)[1:(length(fty)/2)]))
par(mfcol=c(2,1))
plot(seq(-50,50,by=.01),y,type="l",ylab = "Cosine Function")
plot(abs(fty),xlim=c(x-30,x+30),type="l",ylab="Spectral Density in hz")


I know I can compute the frequency by taking the bin value and diving it by the size of the interval(total time of the domain). Since the bins started at 1, when it should be zero, it would thus be
frequency=(BinValue-1)/MaxTime, which does get me the 1/9'th I have in the function above.



I have two quick questions:



First) I am having trouble computing the phase. the density function peaks at 12 (see bottom graph), shouldn't the value then the value of the phase be 2*pi+angle[12]
but I am getting a value of




angle[12]
[1] -2.558724




which puts the phase at 2*pi+angle[12]=3.724462. But that's wrong the phase should be 32 radians. What am I doing wrong?



Second) How do I convert the amplitude which is abs(fty) into the constant in my consine function above, which is 7? abs(fty)[12]=34351.41 , how do I convert this number to 7?. I appreciate your help!










share|cite|improve this question











$endgroup$












  • $begingroup$
    So, by $R$, you mean the statistical analysis package, and not the real numbers?
    $endgroup$
    – Thomas Andrews
    Sep 1 '14 at 17:37










  • $begingroup$
    Yes, the stat program. sorry about the confusion, i'll make the edit. thought this forum would be more apt than stack overflow because it really is more of a math question.
    $endgroup$
    – jessica
    Sep 1 '14 at 17:40














0












0








0


2



$begingroup$


This is a math question with programming application. So I am trying to find 3 things given a certain function in the x domain when transformed into the spectral domain.
1) the Amplitude
2) The Frequency
3) The Phase
In R (statistical software) I have coded the following function:



y=7*cos(2*pi*(seq(-50,50,by=.01)*(1/9))+32)
fty=fft(y,inverse=F)
angle=atan2(Im(fty), Re(fty))
x=which(abs(fty)[1:(length(fty)/2)]==max(abs(fty)[1:(length(fty)/2)]))
par(mfcol=c(2,1))
plot(seq(-50,50,by=.01),y,type="l",ylab = "Cosine Function")
plot(abs(fty),xlim=c(x-30,x+30),type="l",ylab="Spectral Density in hz")


I know I can compute the frequency by taking the bin value and diving it by the size of the interval(total time of the domain). Since the bins started at 1, when it should be zero, it would thus be
frequency=(BinValue-1)/MaxTime, which does get me the 1/9'th I have in the function above.



I have two quick questions:



First) I am having trouble computing the phase. the density function peaks at 12 (see bottom graph), shouldn't the value then the value of the phase be 2*pi+angle[12]
but I am getting a value of




angle[12]
[1] -2.558724




which puts the phase at 2*pi+angle[12]=3.724462. But that's wrong the phase should be 32 radians. What am I doing wrong?



Second) How do I convert the amplitude which is abs(fty) into the constant in my consine function above, which is 7? abs(fty)[12]=34351.41 , how do I convert this number to 7?. I appreciate your help!










share|cite|improve this question











$endgroup$




This is a math question with programming application. So I am trying to find 3 things given a certain function in the x domain when transformed into the spectral domain.
1) the Amplitude
2) The Frequency
3) The Phase
In R (statistical software) I have coded the following function:



y=7*cos(2*pi*(seq(-50,50,by=.01)*(1/9))+32)
fty=fft(y,inverse=F)
angle=atan2(Im(fty), Re(fty))
x=which(abs(fty)[1:(length(fty)/2)]==max(abs(fty)[1:(length(fty)/2)]))
par(mfcol=c(2,1))
plot(seq(-50,50,by=.01),y,type="l",ylab = "Cosine Function")
plot(abs(fty),xlim=c(x-30,x+30),type="l",ylab="Spectral Density in hz")


I know I can compute the frequency by taking the bin value and diving it by the size of the interval(total time of the domain). Since the bins started at 1, when it should be zero, it would thus be
frequency=(BinValue-1)/MaxTime, which does get me the 1/9'th I have in the function above.



I have two quick questions:



First) I am having trouble computing the phase. the density function peaks at 12 (see bottom graph), shouldn't the value then the value of the phase be 2*pi+angle[12]
but I am getting a value of




angle[12]
[1] -2.558724




which puts the phase at 2*pi+angle[12]=3.724462. But that's wrong the phase should be 32 radians. What am I doing wrong?



Second) How do I convert the amplitude which is abs(fty) into the constant in my consine function above, which is 7? abs(fty)[12]=34351.41 , how do I convert this number to 7?. I appreciate your help!







fourier-analysis






share|cite|improve this question















share|cite|improve this question













share|cite|improve this question




share|cite|improve this question








edited Sep 1 '14 at 17:40







jessica

















asked Sep 1 '14 at 17:27









jessicajessica

3831823




3831823












  • $begingroup$
    So, by $R$, you mean the statistical analysis package, and not the real numbers?
    $endgroup$
    – Thomas Andrews
    Sep 1 '14 at 17:37










  • $begingroup$
    Yes, the stat program. sorry about the confusion, i'll make the edit. thought this forum would be more apt than stack overflow because it really is more of a math question.
    $endgroup$
    – jessica
    Sep 1 '14 at 17:40


















  • $begingroup$
    So, by $R$, you mean the statistical analysis package, and not the real numbers?
    $endgroup$
    – Thomas Andrews
    Sep 1 '14 at 17:37










  • $begingroup$
    Yes, the stat program. sorry about the confusion, i'll make the edit. thought this forum would be more apt than stack overflow because it really is more of a math question.
    $endgroup$
    – jessica
    Sep 1 '14 at 17:40
















$begingroup$
So, by $R$, you mean the statistical analysis package, and not the real numbers?
$endgroup$
– Thomas Andrews
Sep 1 '14 at 17:37




$begingroup$
So, by $R$, you mean the statistical analysis package, and not the real numbers?
$endgroup$
– Thomas Andrews
Sep 1 '14 at 17:37












$begingroup$
Yes, the stat program. sorry about the confusion, i'll make the edit. thought this forum would be more apt than stack overflow because it really is more of a math question.
$endgroup$
– jessica
Sep 1 '14 at 17:40




$begingroup$
Yes, the stat program. sorry about the confusion, i'll make the edit. thought this forum would be more apt than stack overflow because it really is more of a math question.
$endgroup$
– jessica
Sep 1 '14 at 17:40










1 Answer
1






active

oldest

votes


















0












$begingroup$

There are several reasons why you don't get what you expect to see:




  • Your definition of the phase is not what you get from the FFT. The phase you get from the FFT is the phase of the first point of the time domain signal, which in your case is
    $$frac{-2 picdot 50}{9}+32=-2.91text{ rad}$$
    The reason why you only get an approximate value is answered in the next point.


  • The frequency you chose does not lie exactly on one of the FFT bins. This is why you get spectral leakage and neither the magnitude nor the phase will be exact estimates of your original signal. Another way to see this is that within the time window there is not an integer number of periods of your signal.


  • The magnitude of the output of the FFT for a sinusoidal input signal is scaled by $N/2$ (where $N$ is the length of the FFT). So your magnitude estimate is $34351cdot 2/10001=6.87$ which is a reasonable value given the spectral leakage mentioned above.



If you want to generate a sinusoidal signal which lies exactly on an FFT bin you can use the following formula:



$$y(n)=Acosleft(2pi nfrac{i}{N}+phiright),quad n=0,1,ldots,N-1;quad i=0,1,ldots N/2$$



where $n$ is the discrete time index, $N$ is the FFT length (assumed even), $i$ is the desired bin index (note that the index starts at $0$), and $A$ and $phi$ are some arbitrary constants. In this case the FFT will give you exact values for $A$ and $phi$. Note that for $phi$ you get of course the principal value $phiin [-pi,pi)$.






share|cite|improve this answer









$endgroup$













    Your Answer





    StackExchange.ifUsing("editor", function () {
    return StackExchange.using("mathjaxEditing", function () {
    StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
    StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
    });
    });
    }, "mathjax-editing");

    StackExchange.ready(function() {
    var channelOptions = {
    tags: "".split(" "),
    id: "69"
    };
    initTagRenderer("".split(" "), "".split(" "), channelOptions);

    StackExchange.using("externalEditor", function() {
    // Have to fire editor after snippets, if snippets enabled
    if (StackExchange.settings.snippets.snippetsEnabled) {
    StackExchange.using("snippets", function() {
    createEditor();
    });
    }
    else {
    createEditor();
    }
    });

    function createEditor() {
    StackExchange.prepareEditor({
    heartbeatType: 'answer',
    autoActivateHeartbeat: false,
    convertImagesToLinks: true,
    noModals: true,
    showLowRepImageUploadWarning: true,
    reputationToPostImages: 10,
    bindNavPrevention: true,
    postfix: "",
    imageUploader: {
    brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
    contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
    allowUrls: true
    },
    noCode: true, onDemand: true,
    discardSelector: ".discard-answer"
    ,immediatelyShowMarkdownHelp:true
    });


    }
    });














    draft saved

    draft discarded


















    StackExchange.ready(
    function () {
    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f916144%2fcomputing-fast-fourier-transform-in-r%23new-answer', 'question_page');
    }
    );

    Post as a guest















    Required, but never shown

























    1 Answer
    1






    active

    oldest

    votes








    1 Answer
    1






    active

    oldest

    votes









    active

    oldest

    votes






    active

    oldest

    votes









    0












    $begingroup$

    There are several reasons why you don't get what you expect to see:




    • Your definition of the phase is not what you get from the FFT. The phase you get from the FFT is the phase of the first point of the time domain signal, which in your case is
      $$frac{-2 picdot 50}{9}+32=-2.91text{ rad}$$
      The reason why you only get an approximate value is answered in the next point.


    • The frequency you chose does not lie exactly on one of the FFT bins. This is why you get spectral leakage and neither the magnitude nor the phase will be exact estimates of your original signal. Another way to see this is that within the time window there is not an integer number of periods of your signal.


    • The magnitude of the output of the FFT for a sinusoidal input signal is scaled by $N/2$ (where $N$ is the length of the FFT). So your magnitude estimate is $34351cdot 2/10001=6.87$ which is a reasonable value given the spectral leakage mentioned above.



    If you want to generate a sinusoidal signal which lies exactly on an FFT bin you can use the following formula:



    $$y(n)=Acosleft(2pi nfrac{i}{N}+phiright),quad n=0,1,ldots,N-1;quad i=0,1,ldots N/2$$



    where $n$ is the discrete time index, $N$ is the FFT length (assumed even), $i$ is the desired bin index (note that the index starts at $0$), and $A$ and $phi$ are some arbitrary constants. In this case the FFT will give you exact values for $A$ and $phi$. Note that for $phi$ you get of course the principal value $phiin [-pi,pi)$.






    share|cite|improve this answer









    $endgroup$


















      0












      $begingroup$

      There are several reasons why you don't get what you expect to see:




      • Your definition of the phase is not what you get from the FFT. The phase you get from the FFT is the phase of the first point of the time domain signal, which in your case is
        $$frac{-2 picdot 50}{9}+32=-2.91text{ rad}$$
        The reason why you only get an approximate value is answered in the next point.


      • The frequency you chose does not lie exactly on one of the FFT bins. This is why you get spectral leakage and neither the magnitude nor the phase will be exact estimates of your original signal. Another way to see this is that within the time window there is not an integer number of periods of your signal.


      • The magnitude of the output of the FFT for a sinusoidal input signal is scaled by $N/2$ (where $N$ is the length of the FFT). So your magnitude estimate is $34351cdot 2/10001=6.87$ which is a reasonable value given the spectral leakage mentioned above.



      If you want to generate a sinusoidal signal which lies exactly on an FFT bin you can use the following formula:



      $$y(n)=Acosleft(2pi nfrac{i}{N}+phiright),quad n=0,1,ldots,N-1;quad i=0,1,ldots N/2$$



      where $n$ is the discrete time index, $N$ is the FFT length (assumed even), $i$ is the desired bin index (note that the index starts at $0$), and $A$ and $phi$ are some arbitrary constants. In this case the FFT will give you exact values for $A$ and $phi$. Note that for $phi$ you get of course the principal value $phiin [-pi,pi)$.






      share|cite|improve this answer









      $endgroup$
















        0












        0








        0





        $begingroup$

        There are several reasons why you don't get what you expect to see:




        • Your definition of the phase is not what you get from the FFT. The phase you get from the FFT is the phase of the first point of the time domain signal, which in your case is
          $$frac{-2 picdot 50}{9}+32=-2.91text{ rad}$$
          The reason why you only get an approximate value is answered in the next point.


        • The frequency you chose does not lie exactly on one of the FFT bins. This is why you get spectral leakage and neither the magnitude nor the phase will be exact estimates of your original signal. Another way to see this is that within the time window there is not an integer number of periods of your signal.


        • The magnitude of the output of the FFT for a sinusoidal input signal is scaled by $N/2$ (where $N$ is the length of the FFT). So your magnitude estimate is $34351cdot 2/10001=6.87$ which is a reasonable value given the spectral leakage mentioned above.



        If you want to generate a sinusoidal signal which lies exactly on an FFT bin you can use the following formula:



        $$y(n)=Acosleft(2pi nfrac{i}{N}+phiright),quad n=0,1,ldots,N-1;quad i=0,1,ldots N/2$$



        where $n$ is the discrete time index, $N$ is the FFT length (assumed even), $i$ is the desired bin index (note that the index starts at $0$), and $A$ and $phi$ are some arbitrary constants. In this case the FFT will give you exact values for $A$ and $phi$. Note that for $phi$ you get of course the principal value $phiin [-pi,pi)$.






        share|cite|improve this answer









        $endgroup$



        There are several reasons why you don't get what you expect to see:




        • Your definition of the phase is not what you get from the FFT. The phase you get from the FFT is the phase of the first point of the time domain signal, which in your case is
          $$frac{-2 picdot 50}{9}+32=-2.91text{ rad}$$
          The reason why you only get an approximate value is answered in the next point.


        • The frequency you chose does not lie exactly on one of the FFT bins. This is why you get spectral leakage and neither the magnitude nor the phase will be exact estimates of your original signal. Another way to see this is that within the time window there is not an integer number of periods of your signal.


        • The magnitude of the output of the FFT for a sinusoidal input signal is scaled by $N/2$ (where $N$ is the length of the FFT). So your magnitude estimate is $34351cdot 2/10001=6.87$ which is a reasonable value given the spectral leakage mentioned above.



        If you want to generate a sinusoidal signal which lies exactly on an FFT bin you can use the following formula:



        $$y(n)=Acosleft(2pi nfrac{i}{N}+phiright),quad n=0,1,ldots,N-1;quad i=0,1,ldots N/2$$



        where $n$ is the discrete time index, $N$ is the FFT length (assumed even), $i$ is the desired bin index (note that the index starts at $0$), and $A$ and $phi$ are some arbitrary constants. In this case the FFT will give you exact values for $A$ and $phi$. Note that for $phi$ you get of course the principal value $phiin [-pi,pi)$.







        share|cite|improve this answer












        share|cite|improve this answer



        share|cite|improve this answer










        answered Sep 2 '14 at 9:03









        Matt L.Matt L.

        8,871822




        8,871822






























            draft saved

            draft discarded




















































            Thanks for contributing an answer to Mathematics Stack Exchange!


            • Please be sure to answer the question. Provide details and share your research!

            But avoid



            • Asking for help, clarification, or responding to other answers.

            • Making statements based on opinion; back them up with references or personal experience.


            Use MathJax to format equations. MathJax reference.


            To learn more, see our tips on writing great answers.




            draft saved


            draft discarded














            StackExchange.ready(
            function () {
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f916144%2fcomputing-fast-fourier-transform-in-r%23new-answer', 'question_page');
            }
            );

            Post as a guest















            Required, but never shown





















































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown

































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown







            Popular posts from this blog

            Aardman Animations

            Are they similar matrix

            “minimization” problem in Euclidean space related to orthonormal basis