Why python np.fft.fft() result and C# Fourier.Forward() result different?

Question

I am trying to make graphs from I,Q data.

Input is an array of complex numbers created by combining I and Q from I Q data in both C# and Python.

• Python: signal.append(complex( I[i],Q[i]))
• C#: signal[i] = new Complex(I[i], Q[i])

However, the results after fft are different.

• Python : X=np.fft.fft(signal)

• C# : Fourier.Forward(signal, FourierOptions.Default);

  **Python result after FFT**
  
  134.99774296753057+123.375408089594j  
  0j    
  -4.1321113197767545e-15-1.1168149738338684e-14j   
  0j    
  1.2906342661267445e-15+1.2220259526518618e-14j
  0j    
  -7.181755190543981e-15-2.5777990853015353e-15j
  0j    
  
  
  **C# result after FFT**
  {(0.603728260168885, 0.55175159848022)}   
  {(5.03413889952152E-15, 6.88380587394808E-13)}    
  {(-1.93271009179521E-12, 1.21853222667336E-12)}   
  {(3.28118559584034E-13, -1.75079741448138E-12)}
  {(-6.73242997306789E-13, -6.27480109322981E-13)}  
  {(-7.67671396112748E-14, -1.17023881533025E-12)}  
  {(1.38179442802571E-12, 3.75998104929882E-13)}    
  {(-3.81679745939109E-13, 8.32395031201847E-13)}
  

Why?

I don't know why there is a 0j in every even position.

What should I do to make the results the same?

I want to create a Frequency domain and Spectrogram graph.

Frequency domain from I Q data Spectrogram from I Q data Complex number array input FFT

Why python and C# FFT result are different?

import numpy as np
import scipy
import matplotlib.pyplot as plt
import math
import scipy.fftpack


samplerate =10000
duration = 0.05
totalsignalduration = 0.5
startFreq = 0
endFreq = 2000
repeat = 10


numSignalSamples = int(samplerate*totalsignalduration)
time =[0]*numSignalSamples
chirpsignal = [0]* (numSignalSamples*repeat)


FreqDelta = (endFreq - startFreq) / (samplerate*duration)

for nRepeat in range(repeat):
    for i in range(numSignalSamples):
        t=i/samplerate
        time.append(t)
        instantFreq = startFreq+(endFreq-startFreq)*t/duration

        if time[i]<=duration:
            instantFreq=startFreq+FreqDelta*i
        else:
            instantFreq=0

        phase=2*math.pi*instantFreq*t

        if time[i]<=duration:
            chirpsignal[i+nRepeat*numSignalSamples] = complex(math.cos(phase), math.sin(phase))
        else:
            chirpsignal[i+nRepeat*numSignalSamples] = complex(0,0)


# print(chirpsignal[0:10])


I = [0]* len(chirpsignal)
Q = [0]* len(chirpsignal)
signal =np.empty_like(chirpsignal, dtype=np.complex128)

for i in range(len(chirpsignal)):
    I[i]= chirpsignal[i].real
    Q[i] =chirpsignal[i].imag
    signal[i]= complex(I[i],Q[i]) 


fftResult = np.fft.fft(signal)

frequencies = []

for i in range(len(signal)):
    frequencies.append(i*samplerate/len(I)/1000)

plt.figure()
plt.plot(frequencies,fftResult)
plt.show()

C# code

private void PlotSpectrum(Complex[] signal, ChirpProperty chirpProperty)
        {
            // Get a reference to the GraphPane
            GraphPane graphPane = SpectrumGraph.GraphPane;
            graphPane.CurveList.Clear();
            double[] frequencies = new double[signal.Length];
            // Create a list of points to plot (sine wave)
            PointPairList pointList = new PointPairList();
            for (int i = 0; i < signal.Length; i++)
            {
                frequencies[i] = i * chirpProperty.SamplingRate / signal.Length;
                pointList.Add(frequencies[i], signal[i].MagnitudeSquared());
            }


            // Add a curve to the graph
            LineItem curve = graphPane.AddCurve(" Spectrum", pointList, System.Drawing.Color.Blue, SymbolType.None);

            // Set axis labels
            graphPane.XAxis.Title.Text = "Freq";
            graphPane.YAxis.Title.Text = "Amplitude";

            // Refresh the graph
            SpectrumGraph.AxisChange();
            SpectrumGraph.Invalidate();
        }



private void btnGenerateChirpIQ_Click(object sender, EventArgs e)
        {

            try
            {
                ChirpSignalGenerator chirpGen = new ChirpSignalGenerator();
                FFTOperation fft = new FFTOperation();
                
                //chirpGen.Generate(ref I, ref Q, chirpProperty);
                chirpGen.Generate(ref srcI, ref srcQ,chirpProperty);


                PlotTime(srcI, srcQ, chirpProperty);
                spectrum = fft.FFT(srcI, srcQ);
                PlotSpectrum( spectrum, chirpProperty );
            }
            catch( Exception ex )
            {
                MessageBox.Show(ex.Message);
            }
        }



class ChirpSignalGenerator
    {
        public StringBuilder sb1 = new StringBuilder();
        public StringBuilder sb2 = new StringBuilder();
        public ChirpSignalGenerator()
        {

        }
        public void Generate( ref double[] I, ref double[] Q, ChirpProperty chirpProperty)
        {
            double sampleRate = chirpProperty.SamplingRate;
            //int sampleRate = 1000; // Sample rate in Hz
            double duration = chirpProperty.Duration; // Duration of the signal in seconds
            double TotalSignalDuration = chirpProperty.TotalSignalDuration;
            double startFrequency = chirpProperty.StartFreq; // Start frequency in Hz
            double endFrequency = chirpProperty.StopFreq; // End frequency in Hz
            int numSignalSamples = (int)(sampleRate * TotalSignalDuration);
            double[] time = new double[numSignalSamples];
            Complex[] chirpSignal = new Complex[numSignalSamples *  chirpProperty.Repeat];
            double FreqDelta = (endFrequency - startFrequency) / (sampleRate * duration);
            for (int nRepeat = 0; nRepeat < chirpProperty.Repeat; nRepeat++)
            {
                for (int i = 0; i < numSignalSamples; i++)
                {
                    double t = i / (double)sampleRate;
                    time[i] = t;
                    double instantFrequency = startFrequency + (endFrequency - startFrequency) * t / duration;
                    instantFrequency = (time[i] <= duration) ? startFrequency + FreqDelta * i : 0;
                    double phase = 2 * Math.PI * instantFrequency * t;
                    chirpSignal[i + nRepeat * numSignalSamples ] = (time[i] <= duration) ? new Complex(Math.Cos(phase), Math.Sin(phase)) : new Complex(0, 0);
                }
            }

            // You now have a complex chirp signal in the chirpSignal array.
            // You can use this array to work with the I and Q components separately.

            // For example, to access the In-phase (I) and Quadrature (Q) components:
            I = new double[chirpSignal.Length];
            Q = new double[chirpSignal.Length];

            for (int i = 0; i < chirpSignal.Length; i++)
            {
                I[i] = chirpSignal[i].Real;
                Q[i] = chirpSignal[i].Imaginary;
            }

            // Now you can use the I and Q components as needed.
        }

     

    }


class FFTOperation
    {
        public FFTOperation()
        {
        }
        public Complex[] FFT( double[] I, double[] Q)
        {
            Complex[] signal = new Complex[I.Length];
            for (int i = 0; i < I.Length; i++)
            {
                signal[i] = new Complex(I[i], Q[i]);
            }


            // Perform FFT
            Fourier.Forward(signal, FourierOptions.Default);
            return signal;

        }
    }

Answer 1

I dug around a bit and I managed to identify your problem. It was a scaling issue, which you can control with FourierOptions.

• The FourierOptions.Default does symmetric scaling
• scipy.fft uses asymmetric scaling

When to use what? Using the symmetric and asymmetric one for real signals is kind of the same, since there are no complex numbers included in the raw signal. However, using symmetric for signals that contain complex values is not correct.

Use FourierOptions.Matlab or if possible FourierOptions.AsymmetricScaling. For me the first option worked and the second threw an error.

Small example in Jupyter interactive C#:

#r "nuget:MathNet.Numerics, 4.15.0"
using System;
using System.Numerics;
using System.Linq;
using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;
// displaying a complex array function
static void DisplayComplexArray(Complex[] array)
{
    foreach (var complexNumber in array)
    {
        Console.WriteLine(complexNumber);
    }
}
// generating a chirp, hopefully no mistakes
static Complex[] GenerateIQChirp(double sampleRate, double duration, double startFrequency, double endFrequency)
    {
    int numSamples = (int)(sampleRate * duration);
    double[] time = Enumerable.Range(0, numSamples).Select(i => i / sampleRate).ToArray();

    Complex[] iqChirp = new Complex[numSamples];

    for (int i = 0; i < numSamples; i++)
    {
        double phase = 2 * Math.PI * (startFrequency * time[i] + 0.5 * (endFrequency - startFrequency) * Math.Pow(time[i], 2));
        iqChirp[i] = new Complex(Math.Cos(phase), Math.Sin(phase));
    }
    return iqChirp;
}
double sampleRate = 20; // fs
double duration = 1.0; // duration
double startFrequency = 1; // start frequency
double endFrequency = 5; // end frequency
Complex[] iqChirp = GenerateIQChirp(sampleRate, duration, startFrequency, endFrequency); // generate chirp
// display signals
Console.WriteLine("Original Signal:"); 
DisplayComplexArray(iqChirp);
Fourier.Forward(iqChirp, FourierOptions.Matlab); // !!!!!!! FourierOptions IMPORTANT!!!!!!!!
Console.WriteLine("FFT of Signal:");
DisplayComplexArray(iqChirp);

The python example:

import numpy as np
from scipy.fft import fft
def displayComplexArray(array):
    for complexNumber in array:
        print(complexNumber)
def generateIQChirp(sampleRate, duration, startFrequency, endFrequency):
    numSamples = int(sampleRate * duration)
    time = np.arange(numSamples) / sampleRate
    iqChirp = np.exp(2j * np.pi * (startFrequency * time + 0.5 * (endFrequency - startFrequency) * time**2))
    return iqChirp
sampleRate = 20  # fs
duration = 1.0    # duration
startFrequency = 1  # start frequency
endFrequency = 5    # end frequency
# Generate IQ chirp signal
iqChirp = generateIQChirp(sampleRate, duration, startFrequency, endFrequency)
print("Original Signal:")
displayComplexArray(iqChirp)
fftResult = fft(iqChirp)
print("FFT of Signal:")
displayComplexArray(fftResult)

Finally, the results from C#:

As you can see, they are the same except for some very small deviations.

Some resources for you to check:

C# FourierOptions

Scipy doc: search for the words "asymmetric spectrum"

Sorry in advance if my C# code is not that great, not really my language...