为什么 python np.fft.fft() 结果和 C# Fourier.Forward() 结果不同？

Question

我正在尝试根据 I、Q 数据制作图表。

输入是通过组合 C# 和 Python 中 IQ 数据中的 I 和 Q 创建的复数数组。

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

然而，fft 后的结果不同。

• Python：X=np.fft.fft(信号)

• 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)}
  

为什么？

我不知道为什么每个偶数位置都有一个0j。

我应该怎样做才能使结果相同？

我想创建频域和频谱图。

IQ 数据的频域 IQ 数据的频谱图 复数数组输入 FFT

为什么Python和C# FFT结果不同？

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#代码

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

我查了一下，终于找到了你的问题。这是一个缩放问题，您可以使用

FourierOptions

来控制。

• FourierOptions.Default

  进行对称缩放

• scipy.fft

  使用非对称缩放

什么时候用什么？对真实信号使用对称和非对称信号是相同的，因为原始信号中不包含复数。然而，对包含复数值的信号使用对称是不正确的。

使用

FourierOptions.Matlab

或如果可能的话使用

FourierOptions.AsymmetricScaling

。对我来说，第一个选项有效，第二个选项抛出错误。

Jupyter 交互式 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);

Python 示例：

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)

最后是C#的结果：

如您所见，除了一些非常小的偏差之外，它们是相同的。

一些资源供您查看：

C# FourierOptions

Scipy 文档：搜索单词“非对称频谱”

提前抱歉，如果我的 C# 代码不是那么好，不是我的语言......