00   /**********************************************************************
01   Copyright (c) 2014 Mark C. Erickson
02  
03   Permission is hereby granted, free of charge, to any person obtaining
04   a copy of this software and associated documentation files (the
05   "Software"), to deal in the Software without restriction, including
06   without limitation the rights to use, copy, modify, merge, publish,
07   distribute, sublicense, and/or sell copies of the Software, and to
08   permit persons to whom the Software is furnished to do so, subject to
09   the following conditions:
10  
11   The above copyright notice and this permission notice shall be included
12   in all copies or substantial portions of the Software.
13  
14   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
15   OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
17   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
18   CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
19   TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
20   SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21   **********************************************************************/
22  

     Disclaimer...

23   desc:EQ, Notch Nth-Order, v1.0.1
24  
25   slider1:2000<20,20000,1>Cutoff [Hz]
26   slider2:2<2,8,2>Order [N]
27   slider3:2<0.5,20,0.1>Q
28  
29   @init
30   M = 1;
31   Q = 2;
32  
33  
34   @slider
35  
36   function find_Dk(index stages)
37   (
38   _k = index;
39   _M = stages;
40   2.0 * sin($pi*(2*_k-1)/(2*_M));
41   );
42  
43   slider1 > 19 ? (
44   slider2 > 1 ? (
45  
46   slider3 > 0.5 ? Q = slider3;
47   Q < ( 2.1 * slider1 / srate ) ? ( Q = 2.1 * slider1 / srate );
48   slider3 = Q;
49   sliderchange(slider3);
50  
51   N = slider2; // order N, 2 < N < 8
52   M = N / 2; // # of cascaded elements M, 1 < M < 4
53   w = 2.0 * $pi * slider1 / srate;
54   dE = 2.0 * tan( w / (2.0 * Q)) / sin(w);
55  
56   k = 1;
57   j = 1;
58   Dk = find_Dk(k, M);
59   Ak = (1 + sqr(dE / 2.0)) / (Dk * dE / 2.0);
60   dk = sqrt(( dE * Dk ) / (Ak + sqrt(sqr(Ak) - 1) ));
61   Bk = Dk * (dE / 2.0) / dk;
62   Wk = Bk + sqrt(sqr(Bk) - 1);
63   theta_1k = 2.0 * atan(tan(w/2.0) / Wk);
64  
65   beta1 = 0.5 * (( 1.0 - dk / 2.0 * sin(theta_1k) ) / ( 1.0 + dk / 2.0 * sin(theta_1k) ));
66   gamma1 = ( 0.5 + beta1 ) * cos(theta_1k);
67   alpha1 = 0.5 * (0.5 + beta1) * ((1 - cos(theta_1k))/(1 - cos(w)));
68  
69   M > 1 ? (
70   k = 1;
71   j = 2;
72   theta_2k = 2.0 * atan( Wk * tan(w/2.0) );
73  
74   beta2 = 0.5 * (( 1.0 - dk / 2.0 * sin(theta_2k) ) / ( 1.0 + dk / 2.0 * sin(theta_2k) ));
75   gamma2 = ( 0.5 + beta2 ) * cos(theta_2k);
76   alpha2 = 0.5 * (0.5 + beta2) * ((1 - cos(theta_2k))/(1 - cos(w)));
77   );
78  
79   M > 2 ? (
80   k = 2;
81   j = 1;
82   Dk = find_Dk(k, M);
83   Ak = (1 + sqr(dE / 2.0)) / (Dk * dE / 2.0);
84   dk = sqrt(( dE * Dk ) / (Ak + sqrt(sqr(Ak) - 1) ));
85   Bk = Dk * (dE / 2.0) / dk;
86   Wk = Bk + sqrt(sqr(Bk) - 1);
87   theta_1k = 2.0 * atan(tan(w/2.0) / Wk);
88  
89   beta3 = 0.5 * (( 1.0 - dk / 2.0 * sin(theta_1k) ) / ( 1.0 + dk / 2.0 * sin(theta_1k) ));
90   gamma3 = ( 0.5 + beta3 ) * cos(theta_1k);
91   alpha3 = 0.5 * (0.5 + beta3) * ((1 - cos(theta_1k))/(1 - cos(w)));
92   );
93  
94   M > 3 ? (
95   k = 2;
96   j = 2;
97   theta_2k = 2.0 * atan( Wk * tan(w/2.0) );
98  
99   beta4 = 0.5 * (( 1.0 - dk / 2.0 * sin(theta_2k) ) / ( 1.0 + dk / 2.0 * sin(theta_2k) ));
100   gamma4 = ( 0.5 + beta4 ) * cos(theta_2k);
101   alpha4 = 0.5 * (0.5 + beta4) * ((1 - cos(theta_2k))/(1 - cos(w)));
102   );
103  
104   );
105   );
106  
107   @sample
108   x_1 = spl0;
109   y_1 = 2.0 * ( alpha1 * ( x_1 - 2.0 * cos(w) * x1_1 + x2_1 ) + gamma1 * y1_1 - beta1 * y2_1 );
110   x2_1 = x1_1;
111   x1_1 = x_1;
112   y2_1 = y1_1;
113   y1_1 = y_1;
114   spl0 = y_1;
115  
116   M > 1 ? (
117   x_2 = spl0;
118   y_2 = 2.0 * ( alpha2 * ( x_2 - 2.0 * cos(w) * x1_2 + x2_2 ) + gamma2 * y1_2 - beta2 * y2_2 );
119   x2_2 = x1_2;
120   x1_2 = x_2;
121   y2_2 = y1_2;
122   y1_2 = y_2;
123   spl0 = y_2;
124   );
125  
126   M > 2 ? (
127   x_3 = spl0;
128   y_3 = 2.0 * ( alpha3 * ( x_3 - 2.0 * cos(w) * x1_3 + x2_3 ) + gamma3 * y1_3 - beta3 * y2_3 );
129   x2_3 = x1_3;
130   x1_3 = x_3;
131   y2_3 = y1_3;
132   y1_3 = y_3;
133   spl0 = y_3;
134   );
135  
136   M > 3 ? (
137   x_4 = spl0;
138   y_4 = 2.0 * ( alpha4 * ( x_4 - 2.0 * cos(w) * x1_4 + x2_4 ) + gamma4 * y1_4 - beta4 * y2_4 );
139   x2_4 = x1_4;
140   x1_4 = x_4;
141   y2_4 = y1_4;
142   y1_4 = y_4;
143   spl0 = y_4;
144   );
145