Forums

Dolph-Chebyshev Window (octave+maxima)

Started by Johann Klammer February 14, 2020
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. 
Here is what i've tried:
in maxima:

> 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). and do:
> 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);
and geta filter kernel(in a csvfile... single line..):
> 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?
Nevermind.. it seems the gamma needs2be integer.

On 02/14/2020 10:07 PM, Johann Klammer wrote:
> 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.
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-]
On 02/23/2020 10:54 PM, Richard (Rick) Lyons wrote:
> 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-] >
from the freq vector in the first post? then it must be octaves ifft()... <http://members.aon.at/~aklamme4/scratch/tf.png>
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-]
On 02/25/2020 12:20 AM, Richard (Rick) Lyons wrote:
> 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-] >
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
Hello Johann. You are very welcome.
Sie k&ouml;nnen mich belohnen, indem Sie mir 5 Kilogramm 
Nurnberger Bratwurst schicken. :-)
On 02/25/2020 09:58 AM, Johann Klammer wrote:
> 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