Dolph-Chebyshev Window (octave+maxima)

[Johann Klammer]

I have followed rick lyons article on this
<https://www.dsprelated.com/showarticle/42.php>
,and using the example values he gives, the output matches.

But using different parameters the filter kenel ends up strangely
wiggly. I must be doing something wrong.

N:32;
gamma:48/20;
M:N+1;
alpha:cosh(acosh(10^gamma)/N);
A(m):=abs(alpha*cos(%pi*m/N));
W(m):=((-1)^m)*chebyshev_t(N,A(m));
makelist(''ev(ratsimp(W(n)),numer),n,0,N-1);

then, I copy the vector over to octave(editing out the newlines).

freq= [251.1886431509343, - 105.9392735931824, 0.4193242363116815, - 0.6394550815653592, 0.01522099347604597, 0.3856457167948537, 0.5997410704897626, 0.7320135982679906, 0.8183989191419395, 0.8771985186325758, 0.9183630436608434, 0.9478210284728955, 0.9668451541730011, 1.000472238713574, 0.8156703834037801, 2.638631937232813, - 11.99999999999998, 2.638631937232824, 0.8156703834037801, 1.000472238713574, 0.966845154173002, 0.9478210284728955, 0.9183630436608434, 0.8771985186325758, 0.8183989191419395, 0.7320135982679906, 0.599741070489752, 0.3856457167948537, 0.01522099347604597, - 0.6394550815653592, 0.4193242363116815, - 105.9392735931824]
td=ifft(freq);
td=real(td);
N=32;
te=[td(1)/2,td(2:N),td(1)/2];
tf=te/max(te);
0.05060350341419111,0.09318854394997424,0.0891740552692859,0.180542307461633,0.194132939061669,0.306667930702
4942,0.334105377095744,0.4643736917713142,0.497180244238745,0.6369967107549066,0.6630579828897345,0.800695027
2395862,0.8068400400763767,0.9292125330996996,0.9046336599899251,1,0.9393020597770186,1,0.9046336599899251,0.
9292125330996994,0.8068400400763767,0.8006950272395861,0.6630579828897345,0.6369967107549066,0.49718024423874
5,0.4643736917713142,0.334105377095744,0.3066679307024942,0.194132939061669,0.180542307461633,0.0891740552692
859,0.09318854394997422,0.05060350341419111

plotting this shows wiggles, unsure why.

what am I doing wrong?

[Johann Klammer]

Nevermind.. it seems the gamma needs2be integer.

[Johann Klammer]

Post by Johann Klammer
Nevermind.. it seems the gamma needs2be integer.

or maybe not.
there seems to also have been an off by one somewhere.
there still was that weird dip in the middle..


i'm using now
M:N+1;
alpha:cosh(acosh(10^gamma)/N);
A(m):=alpha*cos(%pi*m/M); //<<here the M
W(m):=chebyshev_t(N,A(m));
append(makelist(''ev(ratsimp(W(n)),numer),n,0,N/2-1),makelist(''ev(ratsimp(W(n)),numer),n,-N/2,-1));

..and shift the impulse resp manually later.

[Richard (Rick) Lyons]

Hello Johann Klammer. I just now saw your posts.
As far as I know, variable 'gamma' does not have to be an integer.

Looking at your Feb. 13 post, using N = 32 and gamma = 48/20, my MATLAB software computes exactly the same frequency vector as your posted "freq" vector. That is a good thing.

When I execute your code to compute your 'tf' vector (the final window sequence) I obtain the same result as my MATLAB code. As far as I can tell your computed 'tf' sequence is the correct final Chebeshev window sequence.

My computation of your code to compute your 'tf' vector produces:

tf =
0.0627, 0.0740, 0.1135, 0.1629, 0.2223, 0.2910, 0.3677, 0.4507,
0.5375, 0.6252, 0.7106, 0.7904, 0.8611, 0.9197, 0.9636, 0.9908,
1.0000, 0.9908, 0.9636, 0.9197, 0.8611, 0.7904, 0.7106, 0.6252,
0.5375, 0.4507, 0.3677, 0.2910, 0.2223, 0.1629, 0.1135, 0.0740,
0.0627.

When I plot your 'tf' window sequence it looks good to me.

[-Rick-]

[Johann Klammer]

Post by Richard (Rick) Lyons
Hello Johann Klammer. I just now saw your posts.
As far as I know, variable 'gamma' does not have to be an integer.
Looking at your Feb. 13 post, using N = 32 and gamma = 48/20, my MATLAB software computes exactly the same frequency vector as your posted "freq" vector. That is a good thing.
When I execute your code to compute your 'tf' vector (the final window sequence) I obtain the same result as my MATLAB code. As far as I can tell your computed 'tf' sequence is the correct final Chebeshev window sequence.
tf =
0.0627, 0.0740, 0.1135, 0.1629, 0.2223, 0.2910, 0.3677, 0.4507,
0.5375, 0.6252, 0.7106, 0.7904, 0.8611, 0.9197, 0.9636, 0.9908,
1.0000, 0.9908, 0.9636, 0.9197, 0.8611, 0.7904, 0.7106, 0.6252,
0.5375, 0.4507, 0.3677, 0.2910, 0.2223, 0.1629, 0.1135, 0.0740,
0.0627.
When I plot your 'tf' window sequence it looks good to me.
[-Rick-]

from the freq vector in the first post?
then it must be octaves ifft()...

<Loading Image...>

[Richard (Rick) Lyons]

Hi Johann.
I think the problem may be your 'chebyshev_t()' command. My MATLAB 'A(m)' sequence is equal to your 'A(m)' sequence, but my MATLAB 'W(m)' frequency-domain sequence is NOT equal to your 'freq' sequence. My 'W(m)' sequence is:

W(m) =
251.1886 -105.9393, 0.4193 -0.6395, 0.0152, 0.3856, 0.5997,
0.7320, 0.8184, 0.8772, 0.9184, 0.9477, 0.9685, 0.9831, 0.9927,
0.9982, 1.0000, 0.9982, 0.9927, 0.9831, 0.9685, 0.9477, 0.9184,
0.8772, 0.8184, 0.7320, 0.5997, 0.3856, 0.0152 -0.6395,
0.4193, -105.9393.

Perhaps you could use the processing in my blog's Step# 5 to compute your 'freq' sequence, and see what happens.

[-Rick-]

[Johann Klammer]

Post by Richard (Rick) Lyons
Hi Johann.
W(m) =
251.1886 -105.9393, 0.4193 -0.6395, 0.0152, 0.3856, 0.5997,
0.7320, 0.8184, 0.8772, 0.9184, 0.9477, 0.9685, 0.9831, 0.9927,
0.9982, 1.0000, 0.9982, 0.9927, 0.9831, 0.9685, 0.9477, 0.9184,
0.8772, 0.8184, 0.7320, 0.5997, 0.3856, 0.0152 -0.6395,
0.4193, -105.9393.
Perhaps you could use the processing in my blog's Step# 5 to compute your 'freq' sequence, and see what happens.
[-Rick-]

Yes, that fixes it. Thank you.

N=32;
M=N+1;
gamma=48/20;
alpha=cosh(acosh(10^gamma)/N);
m=[0:N-1];
W=((-1).^m).*cheby_t(N,abs(alpha*cos(pi.*m./N)));
td=ifft(W);
td=real(td);
te=[td(1)/2,td(2:N),td(1)/2];
tf=te/max(te);

and...

cheby_t.m
function retval = cheby_t (n,v)
if(abs(v)>1)
retval=cosh(n*acosh(v));
else
retval=cos(n*acos(v));
endif
endfunction

[Richard (Rick) Lyons]

Hello Johann. You are very welcome.
Sie können mich belohnen, indem Sie mir 5 Kilogramm
Nurnberger Bratwurst schicken. :-)

[Johann Klammer]

Post by Johann Klammer
cheby_t.m
function retval = cheby_t (n,v)
if(abs(v)>1)
retval=cosh(n*acosh(v));
else
retval=cos(n*acos(v));
endif
endfunction

I believe this one is more correct:

function retval = chebyshev_t (n,v)
if(v>=1)
retval=cosh(n*acosh(v));
else
if(v<=-1)
retval=((-1) .^ n) * cosh(n*acosh(-v));
else
retval=cos(n*acos(v));
endif
endif
endfunction

[Richard (Rick) Lyons]

Hello Johann.
I can only make your Feb 26, 2020 code work if I multiply your final 'retval' sequence by an n-length sequence of 1, -1, 1, -1, ... . That is, in MATLAB code:

retval = (-1).^(0:n-1).*retval;

[Johann Klammer]

Post by Richard (Rick) Lyons
Hello Johann.
retval = (-1).^(0:n-1).*retval;

Because it's a function file to drop somewhere for octave to load.
It's just the T_N function, not your added phase shift.