What about Magnitude mean after Fourier.Forward?

I code like follower:

 public partial class FormTest : Form
{
    public FormTest()
    {
        this.Icon = Icon.ExtractAssociatedIcon(Application.ExecutablePath);
        InitializeComponent();
    }

    public static double[] GetHarmonic(double[] dChannel, int nOrder, ref double dMagnitude, ref double dInitialPhase)
    {
        int nSamplePoints = dChannel.Length;
        Complex32[] HarmonicComplexArr = new Complex32[nSamplePoints];

        for (int i = 0; i < nSamplePoints; i++)
        {
            HarmonicComplexArr[i] = new Complex32((float)dChannel[i], 0);
        }
        Fourier.Forward(HarmonicComplexArr, FourierOptions.Matlab);//FFT
        //frequency spectrum data
        //double[] ArrFreq = new double[nSamplePoints];
        for (int i = 0; i < nSamplePoints; i++)
        {
            //ArrFreq[i] = HarmonicComplexArr[i].Magnitude;

            if (i == nOrder)//i==0 DC component i==1 First Harmonic,i==2 Second Harmonic
            {
                HarmonicComplexArr[i] *= 2;
                dMagnitude = HarmonicComplexArr[i].Magnitude;
                dInitialPhase = HarmonicComplexArr[i].Phase;
            }
            else
            {
                HarmonicComplexArr[i] = new Complex32(0, 0);
            }
        }

        Fourier.Inverse(HarmonicComplexArr, FourierOptions.Matlab);//逆傅里叶变换

        double[] dHarmonic = new double[nSamplePoints];
        for (int i = 0; i < nSamplePoints; i++)
        {
            dHarmonic[i] = HarmonicComplexArr[i].Real;
        }

        return dHarmonic;
    }

    private void FormTest_Load(object sender, EventArgs e)
    {
        int nSamplePointsPerCycle = 720;
        double[] dChannel = new double[nSamplePointsPerCycle];
        for (int i = 0; i < nSamplePointsPerCycle; i++)
        {
            dChannel[i] = 5 * Math.Cos(2 * Math.PI * i / nSamplePointsPerCycle + 0) + 2 * Math.Cos(4 * Math.PI * i / nSamplePointsPerCycle + 0);
        }

        double dMagnitude = 0;
        double dInitialPhase = 0;
        double[] dFirstHarmonic = GetHarmonic(dChannel, 1, ref dMagnitude, ref dInitialPhase);
        lblFirstMagnitude.Text = dMagnitude.ToString();
        double[] dSecondHarmonic = GetHarmonic(dChannel, 2, ref dMagnitude, ref dInitialPhase);
        lblSecondMagnitude.Text = dMagnitude.ToString();

        chtChannel.Series["OriginalSeries"].Points.Clear();
        chtChannel.Series["FirstHarmonicSeries"].Points.Clear();
        chtChannel.Series["SecondHarmonicSeries"].Points.Clear();
        for (int i = 0; i < nSamplePointsPerCycle; i++)
        {
            chtChannel.Series["OriginalSeries"].Points.AddXY(360 * i / nSamplePointsPerCycle, dChannel[i]);
            chtChannel.Series["FirstHarmonicSeries"].Points.AddXY(360 * i / nSamplePointsPerCycle, dFirstHarmonic[i]);
            chtChannel.Series["SecondHarmonicSeries"].Points.AddXY(360 * i / nSamplePointsPerCycle, dSecondHarmonic[i]);
        }
    }
}

The result of the code show in the fig.1.

  1. The Yellow curve is the first harmonic.The red curve is second harmonic.They are all right because I I multiply the complex by two. Such as “HarmonicComplexArr[i] *= 2;” So the first question is why I need to multiply the complex by two?
  2. You can focus on the two lables at the top of the window,which named ‘First Harmoic Magnitude’ and ‘Second Harmonic Magnitude’. And they are be computed by the code “dMagnitude = HarmonicComplexArr[i].Magnitude;”. I don’t know what dMagnitude mean. And I try to change nSamplePointsPerCycle, the dMagnitude will be changed also. I found 2 * dMagnitude/nSamplePointsPerCycle always equal to Harmoic Curve’s Magnitude. So the second question is What the story about this.

Thanks all for your attention.


Fig.1

Question 1:

Fourier.Forward returns a two-sided spectrum, so you need to add the negative frequencies as well. For nOrder the respective negative frequency is in nSamplesPoint-nOrder. In a real symmetric signal they should have the same magnitude, that’s why multiplying by two works in this case - but only in this case.

Question 2:

I’m not entirely sure what you’re referring to as “Magnitude” here. But in general, when looking at amplitudes in the frequency domain, you need to be aware of their scaling. For some reason you’ve chosen to use Matlab conventions, which do not maintain signal power - i.e. the frequency domain curve has different power than the original signal. The default conventions FourierOptions.Default on the other hand do maintain power, i.e. power of both time and frequency domain signals are the same, for both forward and inverse transforms.

If you’re interested in how the power of the signal is composed by the two harmonics, you can either compute the power of the the resulting curves (sum of squares), or you can indeed skip the inverse transform and just look at the powers at the positive and negative frequency of your order (sum of the two squared magnitudes), which should be the same.

In your example I get the following numbers:

  • Original Signal Power: 10440
  • Harmonic 1 Power: 9000 (86%)
  • Harmonic 1 Frequency domain power sum at frequencies 1 and N-1: 9000
  • Harmonic 2 Power: 1440 (14%)
  • Harmonic 2 Frequency domain power sum at frequencies 2 and N-2: 1440

For this I’ve modified your code along the lines of:

// using FourierOptions.Default
dMagnitude = 0;
for (int i = 0; i < nSamplePoints; i++)
{
	//i==0 DC component i==1 First Harmonic,i==2 Second Harmonic
	if (i == nOrder || i == nSamplePoints - nOrder)
	{
		dMagnitude += HarmonicComplexArr[i].Magnitude * HarmonicComplexArr[i].Magnitude;
	}
	else
	{
		HarmonicComplexArr[i] = Complex32.Zero;
	}
}

Or are you more interested in the time-domain amplitudes?

Does that help?

Thank you very much for your help!
But I also confuse about question 2.
I change the code such as below.

        Fourier.Forward(HarmonicComplexArr, FourierOptions.Matlab);//FFT

        dMagnitude = 0;
        for (int i = 0; i < nSamplePoints; i++)
        {
            if (i == nOrder || i == nSamplePoints - nOrder)//i==0 DC component i==1 First Harmonic,i==2 Second Harmonic
            {
                //HarmonicComplexArr[i] *= 2;
                dMagnitude += HarmonicComplexArr[i].Magnitude;
                dInitialPhase = HarmonicComplexArr[i].Phase;
            }
            else
            {
                HarmonicComplexArr[i] = new Complex32(0, 0);
            }
        }
        dMagnitude /= nSamplePoints;
        dMagnitude *= 2;

        Fourier.Inverse(HarmonicComplexArr, FourierOptions.Matlab);//IFFT

I don’t know why I need add such two codes to get the magnitude in the time-domain.

        dMagnitude /= nSamplePoints;
        dMagnitude *= 2;

Thank you again.

What do you mean with “magnitude in the time domain”? The amplitude?

In the fig.1,
First Harmoic Magnitude is 10, you can see yellow curve.
Second Harmonic Magnitude if 4, you can see red curve.

the original signal is y=5cos(x)+ 2cos(2x);

    for (int i = 0; i < nSamplePointsPerCycle; i++)
    {
        dChannel[i] = 5 * Math.Cos(2 * Math.PI * i / nSamplePointsPerCycle + 0) + 2 * Math.Cos(4 * Math.PI * i / nSamplePointsPerCycle + 0);
    }

You mean the difference between the min and max values of the curves?

Of course.
The magnitude in the time domain is mean the distance of the min and the max values of the curves.
I’ve always understood that. Am I wrong?
Thank you for your patience.

Ah, I see. I’ve never heard of this term being used to describe a curve, but terms and notations differ around the world, no worries. For a periodic signal I usually refer to that as 2 * amplitude, or sometimes “peak to peak amplitude”.

Have a look at the following class:

class Harmonics
{
	int N;
	Complex[] _spectrum;
	
	public Harmonics(double[] signal)
	{
		N = signal.Length;
		_spectrum = Generate.Map(signal, x => new Complex(x, 0.0));
		Fourier.Forward(_spectrum, FourierOptions.AsymmetricScaling);
	}

	public double[] Signal(int order)
	{
		if (order == 0)
		{
			return Generate.Repeat(N, Amplitude(0));
		}
		
		Complex[] s = new Complex[N];
		s[order] = _spectrum[order];
		s[N-order] = _spectrum[N-order];
		Fourier.Inverse(s, FourierOptions.AsymmetricScaling);
		return Generate.Map(s, x => x.Real);
	}

	public double Amplitude(int order)
	{
		if (order == 0)
		{
			return _spectrum[0].Magnitude/N;
		}

		return (_spectrum[order].Magnitude + _spectrum[N - order].Magnitude) / N;
	}

	public double Magnitude(int order)
	{
		if (order == 0)
		{
			return 0;
		}

		return 2*Amplitude(order);
	}

	public double[] FrequencyScale(double sampleRate)
	{
		return Fourier.FrequencyScale(N, sampleRate);
	}

	public double Frequency(int order, double sampleRate)
	{
		return FrequencyScale(sampleRate)[order];
	}
}

I can use this the following way:

int N = 720;
double[] dChannel = new double[N];
for (int i = 0; i < N; i++)
{
        // note: I've added a 3.5 DC offset
	dChannel[i] = 3.5 + 5 * Math.Cos(2 * Math.PI * i / N) + 2 * Math.Cos(4 * Math.PI * i / N + 0);
}

var harmonics = new Harmonics(dChannel);

double a0 = harmonics.Amplitude(0); // 3.5 
double m0 = harmonics.Magnitude(0); // 0.0
double[] s0 = harmonics.Signal(0);

double a1 = harmonics.Amplitude(1); // 5.0
double m1 = harmonics.Magnitude(1); // 10.0
double[] s1 = harmonics.Signal(1);

double a2 = harmonics.Amplitude(2); // 2.0
double m2 = harmonics.Magnitude(2); // 4.0
double[] s2 = harmonics.Signal(2);

double a3 = harmonics.Amplitude(3); // 0.0
double m3 = harmonics.Magnitude(3); // 0.0
double[] s3 = harmonics.Signal(3);

Is this essentially what you’re looking for? Is any of the factor unclear, like the factor of 2 or the division by N?

Excellect job!Will you add the Harmonics Class in Math.Net?
I know that I must * 2 in the first time because of negative frequencies which save in the N - order.
And * 2 in the second time I can’t fully understand the difference of Amplitude and Magnitude.

By the way, how to get the Initial Phase of the Harmonics?
There also positive and negative frequencies.