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
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.
However, the modified code has the same amplitude as the transform to the full range without need to provide t_orig:
(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)
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
.