How to realzie iczt use Matlab?

I'm trying to realize the iczt, but I just find czt function from Matlab. How to realize the iczt, use the czt function or other function?

[John Pickerd]

Here is a link to an 10/11/2019 article where Vladimir Sukhoy and Alexander Stoytchev just took three years of research effort to solve for generalized inverse FFT that includes ICZT. They claim they are the first to solve for this since the invetnion of CZT about 50 years ago. This research was in conjunction with Iowa State University.

https://m.techxplore.com/news/2019-10-year-old-puzzle.html

A more detailed article with the math derivations is located at the following:

https://www.nature.com/articles/s41598-019-50234-9

Answer 1

[John BG]

Hi Wei

try applying twice czt, for instance:

t=[0:.1:10]
f0=.2
x = sin(2*pi*f0*t)
Xcz = czt(x)
x2=-czt(Xcz)
plot(t,101*x,t,x2);grid on

1.- note that the 101 applied on the input signal here, so both signals can be compared, so it's the output that has to be attenuated by length(x) to get back to the signal level.

2.- MATLAB already warns:

Warning: Imaginary parts of complex X and/or Y arguments ignored

So, basically, you get the same signal back, if you ignore the Imaginary part. Rewording, you don't get back the initial signals.

3.- S.Shilling [1] and L.Rabiner R.Schafer [2] point out as CZT limitation that as soon as you take a large amount of samples, (i guess if not using MATLAB vpa, or other tools to handle arbitrary long numbers, loss of precision or) overflow may happen.

If you find this answer of any help solving this question, please click on the thumbs-up vote link, thanks in advance

John

[e-mail]

[1] https://www.google.co.uk/url?sa=t&rct=j&q=&esrc=s&source=web&cd=2&cad=rja&uact=8&ved=0ahUKEwjkkOTeiJzLAhVEEJoKHbvTBxAQFggoMAE&url=https%3A%2F%2Fkrex.k-state.edu%2Fdspace%2Fbitstream%2Fhandle%2F2097%2F7844%2FLD2668R41972S43.pdf&usg=AFQjCNELO9KmErzLpFv0dWKcFlIqSpuUdQ

[2] http://www.ingelec.uns.edu.ar/pds2803/Materiales/Articulos/0086rabi.pdf

[Wei Wang]

 I really appreciate your help.
I try to realize the iczt use your code. I find that if we don't give the w and a ,the matlab will use the following default values:
m = length(x)
w = exp(-j*2*pi/m)
a = 1.
Is it a fft of signal?
And I mean if I have the frequency data like 
Xk = 1 if the k = ω;and Xk = 0 others.Can I get the time domain data use iczt?

[John BG]

Wei

 czt(x,m,w,a)

needs 4 inputs, but

czt(x)

assumes

• m = length(x)
• w = exp(-j*2*pi/m)
• a = 1

I don't know the x you use, i just took a test signal

t=[0:.1:10]
f0=.2
x = sin(2*pi*f0*t)

if you start x with a certain length, first transform

X=czt(x)

and second recover

x2=-czt(czt(x)./length(x))

the imaginary part is not the same, which may be a problem, but the real part still looks like the original signal, just plot as answered

plot(t,101*x,t,x2);grid on

However, if you apply a different m when attempting to recover signal, you may be adding resolution, but you may be adding noise or simply losing data.

Are you attempting to improve frequency resolution in a certain band?

Can you make available a sample of the data you want to process with czt?

Answer 2

[Chris Turnes]

There's not necessarily going to be a straightforward inverse for the Chirp-Z transform unless its parameters define evenly spaced nodes on the unit circle in the complex plane (that is, a = 1 and w = exp(-1j*2*pi/m) for some arbitrary theta and x = length(m). If you do have that, then you can calculate the inverse with

bsxfun(@times, a.^((0:(m-1))'), 1/m*czt(x,m,conj(w),1))

but I'm guessing not since this seems like an odd way to use the CZT.

You can think of czt(x,m,w,a) as applying a matrix A to the input vector x, where the entries of A are A_{i,j} = (a/w^(i-1))^(-(j-1)). This is a Vandermonde matrix, and in general is probably very ill-conditioned unless you happen to have a "nice" problem.

There are formulas for the inverse of a Vandermonde matrix, but this is probably not really going to help you if the matrix is ill-conditioned.

If m and length(x) aren't too large, then you might be willing to explicitly build the matrix and solve using \. Using some random numbers, here's an example:

x = randn(100,1);
a = 0.999;
m = 150;
w = exp(-1j*(2*pi/100 + 2*pi*rand));
% Take the CZT
y = czt(x, m, w, a);
% Try to find the inverse
A = czt(eye(length(x)), m, w, a);
xh = A\y;
norm(x-xh)
ans =
 1.9727e-14

Now, I've made a pretty well-conditioned CZT, so the answer is close. The more ill-conditioned the transform is, the worst off the norm of your answer will be.

If the sizes are large and/or you can't afford (either in time or memory) to construct the matrix, you might instead look to an iterative solver.

Answer 3

[John BG]

Wei

perhaps the following helps, attempt to compare fft and czt:

Fs = 100e6;  % Sampling frequency
fSz = 5000;  % Frame size
A=[1 .1 .01]
f=[5e6 15e6 25e6]
A4=.001;
% f4_low=(f(2)+f(3))/2; f4_high=f(3)+(f(2)+f(3))/2   % too wide sweep
f4_low=f(3)*.8
f4_high=f(3)*1.2

f4=repmat([linspace(f4_low,f4_high,100) linspace(f4_high,f4_low,100)],1,50)   % linear chirp

hsin1 = dsp.SineWave(A(1), f(1), 0, 'SamplesPerFrame', fSz, 'SampleRate', Fs);
hsin2 = dsp.SineWave(A(2), f(2), 0, 'SamplesPerFrame', fSz, 'SampleRate', Fs);
hsin3 = dsp.SineWave(A(3), f(3), 0, 'SamplesPerFrame', fSz, 'SampleRate', Fs);
hsin4 = dsp.SineWave(A4, f4(1), 0, 'SamplesPerFrame', fSz, 'SampleRate', Fs);   % this tone does not have constant frequency

hsb = dsp.SpectrumAnalyzer;
hsb.SampleRate = Fs;
hsb.SpectralAverages = 1;
hsb.PlotAsTwoSidedSpectrum = false;
hsb.RBWSource = 'Auto';
hsb.PowerUnits = 'dBW';
hsb.Position = [749 227 976 721];

y_log=zeros(1,fSz);

for idx = 1:1e4 

%      nph1=pi*randi([0 100],1)/100    % phase noise, when in-phase pulls up the noise floor too high
%      nph2=pi*randi([0 100],1)/100
%      nph3=pi*randi([0 100],1)/100
%      nph4=pi*randi([0 100],1)/100
     y1 = step(dsp.SineWave(A(1), f(1), nph1, 'SamplesPerFrame', fSz, 'SampleRate', Fs));
     y2 = step(dsp.SineWave(A(2), f(2), nph2, 'SamplesPerFrame', fSz, 'SampleRate', Fs));  
     y3 = step(dsp.SineWave(A(3), f(3), nph3, 'SamplesPerFrame', fSz, 'SampleRate', Fs));  
     y4 = step(dsp.SineWave(A4, f4(idx), nph4, 'SamplesPerFrame', fSz, 'SampleRate', Fs)); 
     y5=y1+y2+y3+y4+0.1*randn(fSz,1);       % when noise .0001*randn(fSz,1) everything is easy
     step(hsb,y5);   
     y_log=[y_log;y5'];
end

this is what's displayed when noise .0001*randn(fSz,1)

when you replace .0001*randn(fSz,1) with 0.1*randn(fSz,1);

Now attempt to compare fft to czt:

m=1024
f1=17e6;f2=32e6

Y1abs_maxhold=zeros(1,m)
Y1re_maxhold=zeros(1,m)
Y2abs_maxhold=zeros(1,m)
Y2re_maxhold=zeros(1,m)

for ii=1:1:500

      Y1=fft(y_log(ii,:),m)
      Y1abs=abs(Y1);        % Y1ang=angle(Y1);
      Y1re=real(Y1);          % Y1im=imag(Y1)
      Y1abs_maxhold=[Y1abs_maxhold;max(Y1abs_maxhold(end,:),Y1abs)]
      Y1re_maxhold=[Y1re_maxhold;max(Y1re_maxhold(end,:),Y1abs)]

      w = exp(-j*2*pi*(f2-f1)/(m*Fs));
      a = exp(j*2*pi*f1/Fs);
      Y2 = czt(y_log(ii,:),m,w,a);
      Y2abs=abs(Y2);        % Y2ang=angle(Y2);
      Y2re=real(Y2);          % Y2im=imag(Y2)
      Y2abs_maxhold=[Y2abs_maxhold;max(Y2abs_maxhold(end,:),Y2abs)]
      Y2re_maxhold=[Y2re_maxhold;max(Y2re_maxhold(end,:),Y2abs)]

end

While FFT goes from 0 to whatever max frequency, CZT can be customized to zoom on a narrow frequency range.

I chose f1=17e6 and f2=32e6 around 3rd tone while the fft is flat between 200 and 400 (200 and 400 are not frequencies but fft frequency numerals)

note that a=1 means it does not spiral away or inwards z plane circle radius 1, just zooming on a freq range.

While the FFT catches no 3rd tone

the czt caught something where 3rd (upper) tone expected

Now, back to your original question, is the ICZT possible? so far haven't found any source mentioning any implementation, probably because you need to know arc coefficients in advance, that if not known, have to be swept, and may not hit them at all even if trying to sweep the right arc coefficients range.

Rewording, to reverse the CZT you have to reverse this:

http://www.embedded.com/print/4006427

However, using default values, it is, knowing in advance the length of the sample

y3_rec=real(-czt(mean(Y2abs_maxhold)))
figure(3);plot(y3_rec)
y3_rec(1)=[]                 % 1st sample was ~-1e4
figure(3);plot(y3_rec)
y3_rec=[0 y3_rec]         % keep same length
figure(3);plot(abs(fft(y3_rec)))

insert07

the fft of the back-to-time with -CZT gives 2 tones, which is not clear to tell the frequency, but at least you know there is something there anyway.

my guess: To catch the 4th tone, you may have to sweep different ramps around 3rd tone, and hopefully one of the ramps will show up the 4th tone, but it would take some time to implement such test.

closing joke: If you catch the USS Saratoga leaving Seattle port, don't want to know ;)

If you find this answer of any help solving this question, please click on the thumbs-up vote link,

thanks in advance

John

[e-mail]