0

I'm trying to find the set of fft2 sinusoids that sum up to my original 2d data (entirely real numbers). I found this example, doing the same thing for fft, incredibly helpful (thanks @mark snyder).

I noticed that if I drop the 2nd dimension, then my first column always matches, which tells me I'm not too far off? Finally, why did I have to subtract 1 from input into dimT to align the phase?

Here's what I've tried:

import pandas as pd
import numpy as np
from scipy.fftpack import fftfreq, fft2
import cmath
df=pd.DataFrame(    { 'Ang':[  0.091, -0.141, -0.114], 
                      'Con':[ -0.139, -0.259, -0.573],
                      'EAC':[  0.016, -0.106, -0.044]   } ,
                       columns=['Ang', 'Con', 'EAC']      ,
                       index=[2011,2012,2013]               )

my_fft = fft2(df)
T      , B      =len(df)       , len(df.columns)
freqsT , freqsB = fftfreq(T,1) , fftfreq(B,1)
def ifft2Manually(option):
    recombine      = np.zeros((T,B)) #recombineine   
    for t in list(range(T)):   
        for b in list(range(B)):    
            coef  = abs(my_fft[t][b])/(T*B)
            tfreq , bfreq = freqsT[t] , freqsB[b]
            ps    = cmath.phase(my_fft[t][b])
            dimT =  coef*np.cos(tfreq*2*np.pi*( np.array([df.index]).T - 1  )   + ps  )   
            dimB =  coef*np.cos(bfreq*2*np.pi*( np.array(list(range(B)))    )   + ps  )
            if option==1:
                sinusoid = dimT + dimB
            else: 
                sinusoid=  dimT #exclude dimB here 
            recombine=sinusoid +recombine
            #print('\nt={}\nb={}\ndimT={}\ndimB={}\nsinusoid={}'.format(t,b,dimT,dimB,sinusoid))
    if option==1:
        print('\n nothing matches when both dimensions used') 
    else: 
        print('\n first column matches when ONLY the first dimension (dimT) is used') 
    print('\n original data:\n{} \n\n my attempt to match:\n{}'.format(df, recombine )) 

ifft2Manually(option=1)
ifft2Manually(option=2)

********************** UPDATE: ********************************

This code below now seems to work. My main mistake was that the addition needed to happen inside of the cos() function instead of adding two cos() terms. Hopefully this small working example will help someone else in a similar situation.

import pandas as pd
import numpy as np
from scipy.fftpack import fftfreq, fft2
import cmath

row=2 # arbitrary
col=3 # arbitrary
df=pd.DataFrame(    { 'A':[  0.09, -0.12, -0.11], 
                      'B':[ -0.13,  0.29, -0.57],
                      'C':[  0.01, -0.16, -0.04],
                      'D':[  0.10,  0.65,  0.46]   } ,
                        columns=['A', 'B', 'C','D']   , 
                        index=[2011,2012,2013]               )

print('here is the df:\n', df)
actual = np.array(df)[row][col] 
my_fft2 = fft2(df)
T      , B      =df.shape[0]      , df.shape[1]
mysum= 0 
for t in range(T):   
    for b in range(B): 
        coef  = abs(my_fft2[t][b])/(T*B)
        bfreq , tfreq =  fftfreq(B,1)[b] , fftfreq(T,1)[t]
        ps    = cmath.phase(my_fft2[t][b])
        mysum   += coef*np.cos( 2*np.pi*(tfreq*row + bfreq*col) + ps )   
print('\n replicate value at row={}, col={}'.format( row, col))
print('\n yhat   :', mysum) 
print(' actual :'  , actual )
bri85
  • 111
  • 4
  • 1
    You're using `df.index` where you should be using `np.arange(T)`. Because the index is the year, and that is not used at all when you call `fft2`. `fft2(df)` is the same as `fft2(np.array(df))`. Also, why do you convert the iterator `range` into a list? `for` is designed to work with iterators, if you make it into a list, all you accomplish is two intermediate objects (convert the iterator to a list, then the loop uses an iterator over the list...) – Cris Luengo Aug 02 '23 at 19:02
  • Another question: why would you want to compute the FFT over the rows? I don't know what "Ang", "Con" and "EAC" mean, but I'm sure you can't take these as three samples on a uniform grid of some function. You likely should only compute the FFT along the columns: `np.fft.fft(df, axis=0)`. – Cris Luengo Aug 02 '23 at 19:06
  • Thanks @CrisLuengo! Good points. Honestly, I'm not sure that fft2 is the right approach either, and that's part of the reason I'm trying to deconstruct it a bit so that I can understand it better. Maybe that will the the topic of another post sometime ... – bri85 Aug 02 '23 at 20:12
  • So you'll need to know that `fft2` is the same as applying `fft` to the columns and then to the rows (or the other way around, doesn't matter). It's separable. So you could simplify your inverse code to follow that process, it should be cheaper and easier to write correctly. – Cris Luengo Aug 02 '23 at 20:21
  • Why does `index - 1` work? Because `(2011 - 1) % 3 == 0`. The FFT is periodic with a period of N (which is 3 in your case, you have 3 rows in your example), so instead of `n` in 0..2, you can fill in `n + 3*k`, with any arbitrary integer `k`. – Cris Luengo Aug 02 '23 at 20:27
  • I'm sure it's not a complicated fix, but if someone could help me iron out this toy example, I'd appreciate the help & insight! Cris has already pointed out a couple of my errors. – bri85 Aug 03 '23 at 14:50
  • @CrisLuengo - I'd be very grateful for any insight you can spare! :-) – bri85 Aug 07 '23 at 01:56

0 Answers0