How to realzie iczt use Matlab?
23 views (last 30 days)
Show older comments
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?
1 Comment
John Pickerd
on 11 Oct 2019
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.
A more detailed article with the math derivations is located at the following:
Accepted Answer
John BG
on 29 Feb 2016
Edited: John BG
on 29 Feb 2016
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
2 Comments
John BG
on 29 Feb 2016
Edited: John BG
on 29 Feb 2016
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?
More Answers (2)
Chris Turnes
on 1 Mar 2016
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.
0 Comments
John BG
on 29 Feb 2016
Edited: John BG
on 1 Mar 2016
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173104/image.jpeg)
when you replace .0001*randn(fSz,1) with 0.1*randn(fSz,1);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173105/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173108/image.jpeg)
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173110/image.jpeg)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173111/image.jpeg)
the czt caught something where 3rd (upper) tone expected
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173113/image.jpeg)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173114/image.jpeg)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/173116/image.jpeg)
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
0 Comments
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!