The DFT below can be tuned with an offset to the bins.

This is done at this line, double a = (j+0.5) * sf; where sf is, double sf = 2.0 * Math.PI / n;

j is the bin number index.  To make things simple the I've set the N to 32768 abd the sample rate of the test sine signal to 32768 for a bin width of 1Hz. This aligns the bin number with frequency in Hz.

If the test signal is a 1000.5 Hz adding 0.5 to j in  a = (j+0.5) * sf; centers the 1000th bin to 1000.5Hz

clearing up the error from the the frequency being between bins. This is useful if the signal source is sampled from a fixed frequency analog oscillator. Rather than re-tunning the oscillator frequency to the bin center the bin center can be tuned to the oscillator. Of course this shifts all the bin frequencies.

This works with the DFT but I can't figure out how to do this with an FFT. All attempts have failed.

The FFT is bit reversed type.

private Complex[] Dft(double[] timeSeries)
        {
            UInt32 n = mLengthTotal;
            UInt32 m = mLengthHalf;
            double[] re = new double[m];
            double[] im = new double[m];
            Complex[] result = new Complex[m];
            double sf = 2.0 * Math.PI / n;

            Parallel.For(0, m, (j) =>
            //for (UInt32 j = 0; j < m; j++)
            {
                double a = (j+.5) * sf;

                for (UInt32 k = 0; k < n; k++)
                {
                    re[j] += timeSeries[k] * Math.Cos(a * k) * mDFTScale;
                    im[j] -= timeSeries[k] * Math.Sin(a * k) * mDFTScale;
                }

                result[j] = new Complex(re[j], im[j]);
            });

Why not just mix the input spectrum up or down with e^jw to match the known offset?   It accomplishes the same thing and you don't have to worry about orthogonality issues with the basis functions or making the fast transform work.

Thanks.

I somewhat of a novice with these maths. I'll try to make some sense out of what you said. I'm not sure what you mean by the "input spectrum" or mixing it with the e^jw.

The input spectrum is just the spectrum of the input signal being processed.   It can be shifted up or down in frequency by any arbitrary amount by multiplying the input time sequence by a phasor rotating at the frequency of the desired shift, i.e., to shift by w radians/sec, multiply by e^jwt = cos(wt) + j*sin(wt).   This is often referred to as "mixing" the signal up or down.  Many radio receivers tune the input signal this way.

I was wondering about this. I do understand this in the analog domain.

One would have have to know the frequency offset of the input signal.

What if it's an unknown. The phase of the DFT seems to be fairly linear with the offset from bin center. I was using that to auto correct for the offset. Run the DFT a second time. 

That's a complex pair, (e^jwt = cos(wt) + j*sin(wt). How would I go about multiplying that with real data? 

How would you determine the frequency offset?   Are there assumptions about your input signal you can exploit?   However you detect the offset, you can adjust the shift of the input frequency by adjusting the mix frequency.

The DFT or FFT takes a complex-valued input.   Multiplying either a real-valued input or a complex-valued input by an complex-valued mixer (e^jwt) will result in a complex-valued output that can be processed by a DFT or FFT.

Is necessary to filter the mixed data before the FFT?

Also is there a read on the net about this you can point me to?

One of the beauties of the complex mix (rather than a real-valued mix) is that it translates all of the spectrum equally, so no image-rejection filters are needed.

A web search on terms like "frequency shift complex exponential" should get you started.  Rick Lyons' has a good article here that covers a lot of the details:

https://www.dsprelated.com/showarticle/192.php

Also, look at the frequency shift property of the Fourier Transform, since this is essentially what is being done.

Great. 

Thanks

You might want to create a complex version of you signal by using a Hilbert filter which will then give you a real and imaginary part in the time domain, and then you can shift the signal's frequency as desired, and just keep on the real result for the DFT/FFT.

https://en.wikipedia.org/wiki/Hilbert_transform

I looked at the Wiki for Hilbert transform, and found lots of mathematical terms that you don't need to know for what you are trying to do.

Suffice to say, the Hilbert Filter will give you almost the exact same data as the input, but shifted by 90 degrees in phase over the useful spectrum.  Thus you can treat a delay copy of your input as real and the filtered part as imaginary.

One word of caution.  The real part that is effectively your original signal also has to be delayed by half the length of you Hilbert filter, for the two to line up correctly for 90 degree shift.

Thanks for this.

I though it was going as simple as a just multiplying the input data with a complex pair and doing the shift in one great swoop. I need to be conscious overhead in terms of time and computational demand. Eventually I want to connect all this to a stream of data and run it in win. For now it's just a study. Is it necessary to have the data in IQ form before mixing? I think I just answered that question by asking it.

I believe that you need the complex data input to do the shift correctly.  The Hilbert Filters can be very simple FIR implementation.  The delayed real input is just that simple as to need only a delay.  Just means storage, which incidentally is in the delay line for the Hilbert FIR.

What you are doing is equivalent to using a particular window function.

I have renamed your timeSeries to S to shorten the lines for display purposes.

Here is the core of your code:

  re[j] += S[k] * Math.Cos(a * k) * mDFTScale;

  im[j] -= S[k] * Math.Sin(a * k) * mDFTScale;

By switching to pseudocode for a moment, this can be expressed in complex form:

  result[j] += S[k] * e^(-i * a * k) * mDFTScale;

Substituting in your definition of "a":

   result[j] += S[k] * e^(-i * (j+.5) * sf * k) * mDFTScale;

Separate out the one-half part:

   result[j] += e^(-i * .5 * sf * k) * S[k] * e^(-i * j * sf * k) * mDFTScale;

What you are left with is the standard definition of a DFT with a window function:

   result[j] += w[k] * S[k] * e^(-i * j * sf * k) * mDFTScale;

Where the window function is this:

   w[k] = e^(-i * Pi/n * k)

Therefore, if you want to use a standard FFT, multiply your timeSeries by w[k] and you should get the same results.

Hope this helps,

Ced

Yes. Thanks. It is all helpful. I just have wrap my head around it.

This is a snippet form an FFT routine. The author uses a linked list object rather than an array to hold the real and imaginary data. The x.next statement simply create the next element in the list.

The added code Math.Cos(Math.PI /32768 * k) and Math.Sin(Math.PI /32768 * k) satisfies the frequency shift in the data. All parameter as described in the first post, sample rate etc. Fora 0.5Hz up shift the sample rate in place of n, of Cedron's example, did the trick.

mN is the length of timeSeries;

FFTElement x = mX[0];
UInt32 k = 0;

  for (UInt32 i = 0; i < mN; i++)
  {
    x.re = Math.Cos(Math.PI /32768 * k) * timeSeries[k];
    x.im = Math.Sin(Math.PI /32768 * k) * timeSeries[k];
    x = x.next;
    k++;
  }