

In [1]:

    
import Graphics.Dynamic.Plot.R2
import Data.Complex
import Data.Monoid

We want to express the following series transformation $\operatorname{BFT}: \mathbb{C}^\mathbb{N} \to \mathbb{R}^\mathbb{R} , x\mapsto y$: $$ y(t) = \sum_{k=0}^\infty \left( \Re{x_k} \cdot \cos ((k+\tfrac12)\cdot \pi \cdot t)

                            - \Im{x_k} \cdot \sin ((k+\tfrac12)\cdot \pi \cdot t) \right).



In [2]:

    
plotWindow $ [ plot $ continFnPlot . (.((k+1/2)*pi*))<$>[cos, negate . sin]
             | k <- [0..3] ]
          ++ [forceXRange (-0.5,0.5)]









    





GraphWindowSpecR2{lBound=-0.5, rBound=0.5, bBound=-1.31040238275929, tBound=1.3265344360891878, xResolution=640, yResolution=480}



In [3]:

    
myFourier :: [Complex Double] -> Double->Double
myFourier xs t = sum [ rx * cos ((k+1/2)*pi*t) + ix * sin ((k+1/2)*pi*t)
                     | (k, rx:+ix) <- zip [0..] xs ]

Using $(\cos,i\sin) \varphi = \tfrac12 (e^{i\cdot\varphi} \pm e^{-i\cdot\varphi})$, $$\begin{align} y(t) =& \tfrac12 \cdot \sum_{k=0}^\infty \left( \Re{x_k} \cdot (e^{i\pi\cdot(k+\tfrac12)\cdot t} + e^{-i\pi\cdot(k+\tfrac12)\cdot t})

            - i\Im{x_k} \cdot (e^{i\pi\cdot(k+\tfrac12)\cdot t} - e^{-i\pi\cdot(k+\tfrac12)\cdot t}) \right)

\=& \tfrac12 \cdot \sum_{k=0}^\infty \left( (\Re{x_k} - i\Im{x_i}) \cdot e^{i\pi\cdot(k+\tfrac12)\cdot t}

             + (\Re{x_k} + i\Im{x_i}) \cdot e^{-i\pi\cdot(k+\tfrac12)\cdot t}
             \right)

\=& \sum_{k=0}^\infty \Re\left( xk \cdot (e^{i\pi\cdot(k+\tfrac12)\cdot t})^\ast \right) \=& \Re\left( e^{-i\frac{\pi\cdot t}2} \sum{k=0}^\infty x_k \cdot e^{-i\pi\cdot k\cdot t} \right) \end{align}$$



In [4]:

    
-- myFourier xs t = realPart $ cis (-pi*t/2) * sum [x * cis (-pi*k*t) | (k,x) <- zip [0..] xs]

Discretised over points $-1+\tfrac1n , -1 + \tfrac3n \ldots 1 - \tfrac1{n}$ (midpoints of $n$ subdivisions of the interval $[-1,1]$), i.e. $$\begin{align} t_j = -1 + \tfrac{1+2\cdot j}{n} = \tfrac{1-n}{n} + \tfrac{2\cdot j}{n} \end{align}$$ we get $$\begin{align} y_j =& \Re\left( e^{-\tfrac{i\pi}2\cdot t_j} \sum_{k=0}^{n/2 - 1} x_k \cdot e^{-i\pi\cdot k\cdot \left(\tfrac{1-n}{n} + \tfrac{2\cdot j}{n}\right)} \right) \\=& \Re\left( e^{-\tfrac{i\pi}2\cdot t_j} \sum_{k=0}^{n/2 - 1} x_k \cdot e^{-i\pi\cdot k\cdot\tfrac{1-n}{n}} \cdot e^{-2i\pi\cdot\tfrac{k\cdot j}{n}} \right) \\\equiv& \Re\left( \mu_j \sum_{k=0}^{n - 1} \xi_k \cdot e^{-2i\pi\cdot\tfrac{k\cdot j}{\tilde n}} \right) = \Re\left(\mu_j \cdot \operatorname{DFT}(\xi)_j\right) \end{align}$$

with $$ \xi_k = \begin{cases} x_k \cdot e^{-i\pi\cdot k\cdot\tfrac{1-n}{n}} & \text{for }k<\tfrac{n}2 \\ 0 & \text{else} \end{cases} $$ and $$ \mu_j = e^{-\tfrac{i\pi}2\cdot t_j}. $$ Thus, the transform can be implemented using FFT standards:



In [5]:

    
import Numeric.FFT.Vector.Invertible
import qualified Data.Vector.Storable as SArr
import qualified Data.Vector.Generic as Arr

myFFourier :: [Complex Double] -> [(Double,Double)]
myFFourier xs = [ (t, realPart $ μ*υ)
                | (j,υ) <- zip [0..] $ SArr.toList (run dft ξs)
                , let μ = cis $ -pi * t/2
                      t = ((1-n)/n + 2*j/n)]
 where ξs = SArr.fromList $
               [x * cis (-pi*k*(1-n)/n) | (k,x) <- zip [0..] xs]
             ++ (const 0<$>xs)
       n = fromIntegral $ length xs * 2



In [6]:

    
someCoefs :: [Complex Double]
someCoefs = [0.4:+0.0, 0.6:+1.0, (-0.7):+(-0.6), 0.6:+(-0.9), 0.0:+1.0, 0.1:+0.1, (-0.6):+(-0.2), 0.7:+(-1.0)
            , 0, 0, 0, 0 ]



In [7]:

    
plotWindow [ continFnPlot $ myFourier someCoefs
           , lineSegPlot $ myFFourier someCoefs ]









    





GraphWindowSpecR2{lBound=-1.277777777777778, rBound=1.277777777777778, bBound=-5.318355746967219, tBound=5.851583402618405, xResolution=640, yResolution=480}

To perform the inverse transform, we first need to reconstruct the discarded imaginary parts of $\mu_j \cdot \operatorname{DFT}(\xi)_j$, since that DFT actually operates on $\mathbb{C}^{n}$ rather than $\mathbb{R}^n$. This is possible because only $$ \partial_t \operatorname{BFT}(\rho)(t) = \operatorname{BFT}((i\pi\cdot (k+\tfrac12)\cdot \rho_k)_{(k)})(t) $$ while for real $\iota$ $$ \partial_t \operatorname{BFT}(i\cdot \iota)(t) = -\operatorname{BFT}((\pi\cdot (k+\tfrac12)\cdot \iota_k)_{(k)})(t). $$ We can sort out the real/imaginary parts since they correspond to even/odd components of the signal: for $x = \rho + i\cdot \iota$, $$\begin{align} \operatorname{BFT}(\rho)(t) = \tfrac12 (\operatorname{BFT}(x)(t) + \operatorname{BFT}(x)(-t) \\ \operatorname{BFT}(i\cdot\iota)(t) = \tfrac12 (\operatorname{BFT}(x)(t) - \operatorname{BFT}(x)(-t). \end{align}$$ Furthermore, it is $$\begin{align} \operatorname{BFT}(\rho)(t+1) = \tfrac12 (\operatorname{BFT}(x)(t) + \operatorname{BFT}(x)(-t) \\ \operatorname{BFT}(i\cdot\iota)(t) = \tfrac12 (\operatorname{BFT}(x)(t) - \operatorname{BFT}(x)(-t). \end{align}$$ The imaginary part of the (linear!) DFT must correspond to the real part of the complex-rotated coefficients;



In [8]:

    
myInvFFourier :: Int -> (Double -> Double) -> [Complex Double]
myInvFFourier n' f = [ 2 * cis (pi*k*(1-n)/n) * ξ
                     | (k,ξ) <- zip [0..n/2-1] $ SArr.toList ξs ]
 where ts = SArr.generate n' $ \j' -> (1-n)/n + 2*fromIntegral j'/n
       ys = SArr.map f ts
       μinvs = SArr.map (\t -> cis $ pi*t/2) ts
       ξs = run idft μinvys
       μinvys = SArr.zipWith (*) μinvs $ SArr.map (:+0) ys
       n = fromIntegral n'



In [9]:

    
f x = sin (exp $ 3*x) - x^2
--let ((μinvys, μinvðt'ys),_) = myInvFFourier 8 f
-- in plotWindow [ legendName (pn++"("++vn++")") . lineSegPlot . zip [0..] $ p<$>vs
  --             | (vn,vs) <- [("μ⁻¹·y",μinvys), ("μ⁻¹·∂y/∂t",(2/pi*)<$>μinvðt'ys)]
    --           , (pn,p) <- [("Re",realPart), ("Im",imagPart)] ]
plotWindow $ continFnPlot f
           : [ legendName ("n="++show n) . lineSegPlot . myFFourier $ myInvFFourier n f
             | n<-[4,8,16,32] ]









    





GraphWindowSpecR2{lBound=-1.291666666666667, rBound=1.291666666666667, bBound=-2.132912155457216, tBound=1.4220469953686115, xResolution=640, yResolution=480}



In [ ]: