24 views (last 30 days)

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

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?

Sign in to comment.

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.

Sign in to comment.

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)

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:

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

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/270528-how-to-realzie-iczt-use-matlab#comment_755354

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/270528-how-to-realzie-iczt-use-matlab#comment_755354

Sign in to comment.