Fractional delay using time-shift property of FFT. Why not work?

[pfeatherstone]

I'm experimenting with applying fractional delay using the time-shift property of the FFT. This kind of thing is done inside the DLLs in OFDM demods, particularly 5G.
I wanted to see if this could be done generally rather than convolving with a FIR filter (with appropriate delay) or using a polyphase filter bank.

Here is some code:

#include <vector>
#include <complex>
#include <liquid/liquid.h>
#include "matplotlibcpp.h"

using namespace std;
namespace plt = matplotlibcpp;

/*! HELPERS !*/
vector<float> real(const vector<complex<float>>& z)
{
    vector<float> r(z.size());
    std::transform(z.begin(), z.end(), r.begin(), [](const complex<float>& x) {return x.real();});
    return r;
}

vector<float> imag(const vector<complex<float>>& z)
{
    vector<float> r(z.size());
    std::transform(z.begin(), z.end(), r.begin(), [](const complex<float>& x) {return x.imag();});
    return r;
}
/*! HELPERS !*/

vector<complex<float>> make_example_signal()
{
    int          ftype   = LIQUID_FIRFILT_ARKAISER;
    unsigned int k       =     2;
    unsigned int m       =     9;
    float        beta    = 0.30f;
    int          ms      = LIQUID_MODEM_BPSK;
    symstreamcf  gen     = symstreamcf_create_linear(ftype,k,m,beta,ms);

    vector<complex<float>> iq(1024);
    symstreamcf_write_samples(gen, iq.data(), iq.size());
    symstreamcf_destroy(gen);

    return iq;
}

vector<float> spectrum(const vector<complex<float>>& iq, int nfft)
{
    spgramcf obj = spgramcf_create_default(nfft);
    vector<float> spec(nfft);
    spgramcf_write(obj, (liquid_float_complex*)iq.data(), iq.size());
    spgramcf_get_psd(obj, spec.data());
    spgramcf_destroy(obj);
    return spec;
}

vector<complex<float>> delay1(const vector<complex<float>>& iq, float dt)
{
    const float fc = 0.45;
    vector<float> h((int(ceil(40.0 / (2*fc))) / 2) * 2 + 1); //good length is about ceil(4/rolloff) and a good rolloff is ~ bw/10
    liquid_firdes_kaiser(h.size(), fc, 60.0, dt, h.data());

    vector<complex<float>> out(iq.size() + h.size() - 1);
    firfilt_crcf q = firfilt_crcf_create(h.data(),h.size());
    firfilt_crcf_execute_block(q, (liquid_float_complex*)iq.data(), iq.size(), out.data());
    for (size_t i = iq.size() ; i < out.size() ; i++)
    {
        firfilt_crcf_push(q, 0);
        firfilt_crcf_execute(q, &out[i]);
    }
    firfilt_crcf_destroy(q);

    //gonna remove the filter delay
    std::rotate(out.begin(), out.begin() + (h.size() / 2), out.end());
    out.resize(iq.size());
    return out;
}

vector<complex<float>> delay2(const vector<complex<float>>& iq, float dt)
{
    vector<complex<float>> out(iq.size());
    fft_run(iq.size(), (liquid_float_complex*)iq.data(), out.data(), LIQUID_FFT_FORWARD, 0);
    for (size_t i = 0 ; i < out.size() ; i++)
        out[i] *= std::polar<float>(1, (2*M_PI*dt*i)/out.size());
    fft_run(out.size(), out.data(), out.data(), LIQUID_FFT_BACKWARD, 0);
    const float scale = 1.0 / out.size();
    std::transform(out.begin(), out.end(), out.begin(), [scale](const complex<float>& z) {return scale*z;});
    return out;
}

int main()
{
    float delay = 1.0;

    auto iq1 = make_example_signal();
    auto iq2 = delay1(iq1, delay);
    auto iq3 = delay2(iq1, delay);

    auto spec1 = spectrum(iq1, 512);
    auto spec2 = spectrum(iq2, 512);
    auto spec3 = spectrum(iq3, 512);

    plt::figure();
    plt::named_plot("ori", spec1);
    plt::named_plot("delay - firfilt", spec2);
    plt::named_plot("delay - IFFT(shift(FFT))", spec3);
    plt::legend();

    plt::figure();
    plt::named_plot("ori",real(iq1));
    plt::named_plot("delay - firfilt",real(iq2));
    plt::named_plot("delay - IFFT(shift(FFT))",real(iq3));
    plt::legend();

    plt::figure();
    plt::named_plot("ori",imag(iq1));
    plt::named_plot("delay - firfilt",imag(iq2));
    plt::named_plot("delay - IFFT(shift(FFT))",imag(iq3));
    plt::legend();

    plt::show();

    return 0;
}

Now if I set delay = 1 as above i get the following graphs:

This all looks correct!

Now if I set delay = 0.5, you get the following:

What happened? Why is this gibberish?
Does the time shift property of the FFT only work for integral offsets?

[pfeatherstone]

@jgaeddert What do you think? Has this got to do with the fact the FFT assumes the input is an infinitely repeating sequence with period equal to the FFT size? And for whatever reason, you don't get any nasty effects when the time delay is integer?

Also, it looks like delay2() is essentially doing an FFT convolution with some filter. Indeed

for (size_t i = 0 ; i < out.size() ; i++)
        out[i] *= std::polar<float>(1, (2*M_PI*dt*i)/out.size());

is multiplying the FFT of the input with impulse response H such that H[i] = std::polar<float>(1, (2*M_PI*dt*i)/out.size())

So it's like it's doing an FFT convolution with a filter h such that FFT(h) = H.
Is h a dirac delta with an appropriate shift?
So is it the case that h is only well behaved when dt is an integer?

[jgaeddert]

So, the DFT (and by extension the FFT and IFFT) are cyclic operations. 5G is based on OFDM, and its cyclic prefix allows for a cyclic delay operation such as adjusting the phase slope across frequency from the FFT. But when you take its inverse, you are cyclicly extending its contents. This is fine for OFDM with a cyclic prefix, but with non-cyclic signals (e.g. a single carrier system), you need to account for (and remove) this extension.

[jgaeddert]

To be fair, I didn't scrub through all of your code. But fractional sample offsets are perfectly acceptable with FFT-based time shifting. Looking at your function delay2(), though, I'm not seeing what would be causing the issue unless there's something odd happening with the std::polar() function.

[pfeatherstone]

I doubt std::polar() has any problems, that's just a standard library function. So are you saying FFT-based fractional time shifting ONLY works in OFDM precisely because there is a cyclic prefix and therefore you have a repeating signal? And do you agree that because i don't have a repeating signal, just a BPSK, it is expected that delay2() won't work?

[jgaeddert]

You certainly can use an FFT for adjusting time but you need to zero-pad the FFT and then compensate for the overlap between FFT blocks. Kinda hard to explain just typing on my phone, but if you take a look at fftfilt that might help!

[fsheikh]

The way I understood this is becomes for delay 0.5. The second expression suggest that sequence shift property assumes a DFT periodic in 2*N. With this in mind, I modified the function as below, this is not generic enough and will work for only delay < 0.5, for general case we need to separate the fractional delay in numerator and denominator parts to extend the DFT and then apply the numerator as the integral delay part.

vector<complex<float>> delay2(const vector<complex<float>>& iq, float dt)
{

    size_t period_correction = 1;
    if (dt > 0.0f and dt < 1.0f) {
        // fractional rate desired
        period_correction = static_cast<size_t>(std::round(1/dt));
    }
    vector<complex<float>> out(iq.size() * period_correction, complex<float>(0,0));
    fft_run(iq.size() * period_correction, (liquid_float_complex*)iq.data(), out.data(), LIQUID_FFT_FORWARD, 0);
    for (size_t i = 0 ; i < out.size() ; i++)
        out[i] *= std::polar<float>(1, (2*M_PI*dt*i)/iq.size());
    fft_run(out.size(), out.data(), out.data(), LIQUID_FFT_BACKWARD, 0);
    const float scale = 1.0 / out.size();
    std::transform(out.begin(), out.end(), out.begin(), [scale](const complex<float>& z) {return scale*z;});

    return out;
}