garrettj403/CZT

Zoom-DFT

LarsHuebner-LHNW opened this issue · 1 comments

Hi again @garrettj403

I am not sure if there was a special reason to implement time2freq and freq2time the way you did, therefore I did not create a pull request yet.
I have some modified code which takes care of the proper phase and scaling in time2freq and freq2time, and drops the requirement of providing an original frequency or time. Moreover I think, there is something wrong with the scaling in freq2time if not using the original time, but a zoomed in part.


I put it into a notebook to test everything. Had to zip it, because github does not allow ipynb files:
time2freq_new-zoom-dft.ipynb.zip

Here is a quick overview of its contents:


For time2freq the new code shows the same result as your original code, without the need to compute the f_orig
time2freq_original
time2freq_modified


Here you can see what I mean by the incorrect scaling in freq2time. From a quick test I would guess, that there is a factor of len(t_zoom)/len(t) missing. I could be wrong here, but I looked at the example notebook "example/time-to-freq-to-time-domain-conversion.ipynb" and there the recover of the the original function is only working, because there the time->freq conversion has no f_orig and therefor not the proper scaling, so it cancels the improper scaling in freq->time. Don't know if this is the desired behavior.
freq2time_original

However, the modified code has the same amplitude as the transform to the full range without need to provide t_orig:
freq2time_modified

(and of course units should be ms :) )


def time2freq_new(t, x, f=None):
    """Convert signal from time-domain to frequency-domain.

    Args:
        t (np.ndarray): time
        x (np.ndarray): time-domain signal
        f (np.ndarray): frequency for output signal

    Returns:
        np.ndarray: frequency-domain signal

    """

    # Input time array
    dt = t[1] - t[0]  # time step

    # Output frequency array
    if f is None:
        Nf = len(t)  # number of time points
        # proposed change:
        # f = np.fft.fftshift(np.fft.fftfreq(Nf, dt)) # fftshift -> frequencies in normal ordering
        # for now stay compatible with original code:
        f = np.linspace(-1/(2*dt),1/(2*dt),Nf)
    else:
        Nf = len(f)
        f = np.copy(f)
    df = f[1]-f[0]

    # Step
    W = np.exp(-2j * np.pi * df * dt)

    # Starting point
    A = np.exp(2j * np.pi * f[0] * dt)

    # phase
    phase = np.exp(-2j*np.pi * t[0] * f)

    # Frequency-domain transform
    freq_data = czt.czt(x, Nf, W, A)*phase

    return f, freq_data

def freq2time_new(f, X, t=None):
    """Convert signal from frequency-domain to time-domain.

    Args:
        f (np.ndarray): frequency
        X (np.ndarray): frequency-domain signal
        t (np.ndarray): time for output signal
    Returns:
        np.ndarray: time-domain signal

    """

    # Input frequency array
    df = f[1] - f[0]  # time step

    # Output time array
    if t is None:
        Nt = len(f)  # number of time points
        # proposed change:
        # t = np.fft.fftshift(np.fft.fftfreq(Nt, df)) # fftshift -> time in normal ordering
        # for now stay compatible with original code:
        t = np.linspace(0,1/df,Nt)
    else:
        Nt = len(t)
        t = np.copy(t)
    dt = t[1]-t[0]

    # Step
    W = np.exp(2j * np.pi * dt * df)

    # Starting point
    A = np.exp(-2j * np.pi * t[0] * df)

    # phase
    phase = np.exp(2j*np.pi * f[0] * t)

    # Frequency-domain transform
    # for dft and idft I will use czt, since A and W are chosen properly.
    time_data = czt.czt(X, Nt, W, A)*phase/len(f)

    return t, time_data

The following tests did not show an assertion error, so I think the new functions should work as a replacement for the oririnal ones.

np.testing.assert_almost_equal(
    czt.freq2time(*czt.time2freq(t,sig)),
    freq2time_new(*time2freq_new(t,sig)),
    decimal=12)

np.testing.assert_almost_equal(
    czt.time2freq(t,sig,freq_zoom,f_orig=freq),
    time2freq_new(t,sig,freq_zoom),
    decimal=12)

Hi @LarsHuebner-LHNW

Thank you for pointing this out. I had a bug in the original implementation of time2freq and freq2time. The t_orig and f_orig arguments were a quick and dirty way to fix the bug.

I made a few changes to time2freq and freq2time, based on your recommendations above. They should be working now. I also added a new test to test the phase correction in freq2time.