1%Spencer Jackson 2%a bunch of stuff for a spring reverb 31; 4 5%Time version 6 7%generate an all-pass fractional delay using the Thiran approximation This interpolates between samples to provide a non-integer delay. 8%sys = thiran(tau, Ts) 9%inputs: 10%tau - time of delay 11%Ts - sampling time 12% 13%Outputs: 14%sys - transfer function of the filter 15% 16%Example: 17% H = thiran(4.4,.5); 18% generates 9th order fractional delay filter delaying 4.4 seconds 19% 20% H = thiran(8,.01); 21% generate pure delay of 8 seconds 22%See Also: thiran2 23function sys = thiran(tau, Ts) 24 D = tau/Ts; 25 N = ceil(D); 26 if(N==D) 27 sys = tf([1],[1 zeros(1,N)],Ts); 28 else 29 30 a = (-1).^(1:N).*bincoeff(N,1:N).*prod((D-N+[0:N]'*ones(1,N))./(D-N+(1:N)+[0:N]')); 31 sys = tf([a(N:-1:1) 1],[1 a],Ts); 32 end 33end 34 35%the code below implements the thiran approximation: 36%a = ones(1,N); 37%for k=1:N 38% for i=0:N 39% a(k) = a(k)*(D-N+i)/(D-N+k+i); 40% end 41% a(k) = (-1)^k*bincoeff(N,k)*a(k); 42%end 43 44%sample delay version 45 46%generate an all-pass fractional delay using the Thiran approximation. This interpolates between samples to provide a non-integer delay. 47%sys = thiran(tau, Ts) 48%inputs: 49%D - delay in samples 50%Ts - sampling time 51% 52%Outputs: 53%sys - transfer function of the filter 54% 55%Example: 56% H = thiran2(4.4, 8000); 57% generate 5th order fractional delay filter delaying 4.4 samples 58% 59% H = thiran2(8, 100); 60% generate pure time delay of 8 samples (.08 seconds) 61%See Also: thiran 62function sys = thiran2(D, Ts) 63 N = ceil(D); 64 if (N==D) 65 sys = tf([1],[1 zeros(1,N)],Ts); 66 else 67 a = (-1).^(1:N).*bincoeff(N,1:N).*prod((D-N+[0:N]'*ones(1,N))./(D-N+(1:N)+[0:N]')); 68 sys = tf([a(N:-1:1) 1],[1 a],Ts); 69 end 70end 71 72Fs = 44100/4;%48000; 44100 73Ts = 1/Fs; 74Nyquist = Fs/2; 75 76 77%first lets work on the low frequency chirp section 78%allpass 79Fc = 2800;%4300;% 3700, 2800; 80k = 3.56;%Fs/Fc;%stretch factor this affects Fc. %k=10=4.4khz ,11=4k, 12=3.6k,17=2.6k, 18=2.4k, 9=5k 81Fc = Fs/k/4; 82d = k-floor(k);%larger d makes more stable 83ki = floor(k); 84%for decimal part of k, d E(0,1) a2 = (1-d)/(d+1) 85%then allpass (single stage) is 86%H = (a1 + a1a2z^-1 + a2z^-k + z^-k-1 )/ (1 + a2z^-1 + a1a2z^-k + a1z^-k-1) 87 88al = -.2;%allpass filter coeff (negative due to changes), this stretches the delay curve 89M = 34;% Number of allpass stages (50?) 23 stable limit for al = -.6 90 %k=17;al=-.6;M=23;%2700hz spring 91 %k=12;al=-.36;M=35;3600hz 92 %k=10;al=-.3;M=40;4300hz 93 %k=10;al=-.7;M=19; 94 %k=10;al=-.5;M=30; 95%Hap = (tf([al zeros(1,k-1) 1],[1 zeros(1,k-1) al],Ts))^M;%old way 96a2 = 0; 97if(d) 98 a2 = (1-d)/(d+1); 99end 100Hap = (tf([al al*a2 zeros(1,ki-2) a2 1],[1 a2 zeros(1,ki-2) al*a2 al],Ts))^M; 101 102a = [al al*a2 zeros(1,ki-2) a2 1]; 103b = [1 a2 zeros(1,ki-2) al*a2 al]; 104 105%break 106 107 %so for 3 springs (k al M Fcfreq) 108 %1 17 -.6 23 2594hz 109 %2 12 -.36 35 3675hz 110 %3 10 -.3 40 4410hz 111 %4 10 -.7 19 4410hz 112 %5 10 -.5 30 4410hz 113%new design decimate x4 114%4300hz k=2.564; al=-.7; M=16; D=63ms 115%4300hz k=2.564; al=-.5; M=24; D=88ms 116%4300hz k=2.564; al=-.3; M=33; D=56ms 117%3700hz k=2.9797;al=-3.6;M=42; D=56.2ms 118%2800hz k=3.9375;al=-.6; M=22; D=56.15ms 119 120 121%stretching factor (k) is inversely proportional to Fc 122%Fc = Fs/k;%"cutoff" for a spring, where the maximum dispersion occurs 123if(1) 124 [y t] = impulse(Hap,.06); 125 specgram(y,64,Fs,hanning(64),64-8) 126 %specgram(y,256,Fs,hanning(256),256-32) 127 break 128 %pause; 129end 130break 131 132if(0) 133 k=2.56;al=-.7;M=16;%this way works for any rate 134 %Fc = 4300;%4300;% 3700, 2800; 135 %k = Fs/Fc; 136 d = k-floor(k);%larger d makes more stable 137 ki = floor(k); 138 a2 = 0; 139 if(d) 140 a2 = (1-d)/(d+1); 141 end 142 Hap1 = (tf([al al*a2 zeros(1,ki-2) a2 1],[1 a2 zeros(1,ki-2) al*a2 al],Ts))^M; 143 Hap1 = Hap1*tf(1,[1 zeros(1,ceil(.007*Fs))],Ts); 144 145 k=2.56;al=-.5;M=24; 146 %Fc = 4300;%4300;% 3700, 2800; 147 %k = Fs/Fc; 148 d = k-floor(k);%larger d makes more stable 149 ki = floor(k); 150 a2 = 0; 151 if(d) 152 a2 = (1-d)/(d+1); 153 end 154 Hap2 = (tf([al al*a2 zeros(1,ki-2) a2 1],[1 a2 zeros(1,ki-2) al*a2 al],Ts))^M; 155 Hap2 = Hap2*tf(1,[1 zeros(1,ceil(.032*Fs))],Ts); 156 157 k=2.56;al=-.3;M=33; 158 %Fc = 4300;%4300;% 3700, 2800; 159 %k = Fs/Fc; 160 d = k-floor(k);%larger d makes more stable 161 ki = floor(k); 162 a2 = 0; 163 if(d) 164 a2 = (1-d)/(d+1); 165 end 166 Hap3 = (tf([al al*a2 zeros(1,ki-2) a2 1],[1 a2 zeros(1,ki-2) al*a2 al],Ts))^M; 167 168 k=2.9;al=-.36;M=40; 169 %Fc = 3700;%4300;% 3700, 2800; 170 %k = Fs/Fc; 171 %M = 16;%48k 172 d = k-floor(k);%larger d makes more stable 173 ki = floor(k); 174 a2 = 0; 175 if(d) 176 a2 = (1-d)/(d+1); 177 end 178 Hap4 = (tf([al al*a2 zeros(1,ki-2) a2 1],[1 a2 zeros(1,ki-2) al*a2 al],Ts))^M; 179 Hap4 = Hap4*tf(1,[1 zeros(1,ceil(.0002*Fs))],Ts); 180 181 k=3.9;al=-.6;M=22; 182 %Fc = 2800;%4300;% 3700, 2800; 183 %k = Fs/Fc; 184 %M = 13;%48k 185 %d = k-floor(k);%larger d makes more stable 186 ki = floor(k); 187 a2 = 0; 188 if(d) 189 a2 = (1-d)/(d+1); 190 end 191 Hap5 = (tf([al al*a2 zeros(1,ki-2) a2 1],[1 a2 zeros(1,ki-2) al*a2 al],Ts))^M; 192 Hap5 = Hap5*tf(1,[1 zeros(1,ceil(.00015*Fs))],Ts); 193 194 [y1 t] = impulse(Hap1,.1); 195 [y2 t] = impulse(Hap2,.1); 196 [y3 t] = impulse(Hap3,.1); 197 [y4 t] = impulse(Hap4,.1); 198 [y5 t] = impulse(Hap5,.1); 199 yt = y1+y2+y3+y4+y5; 200 201 specgram(yt,64,Fs,hanning(64),64-8) 202 %specgram(yt,256,Fs,hanning(256),256-32) 203 204 tmp = Hap1.num{1}; 205 save('Spring1Coeff.txt','tmp','-ascii'); 206 tmp = Hap2.num{1}; 207 save('Spring2Coeff.txt','tmp','-ascii'); 208 tmp = Hap3.num{1}; 209 save('Spring3Coeff.txt','tmp','-ascii'); 210 tmp = Hap4.num{1}; 211 save('Spring4Coeff.txt','tmp','-ascii'); 212 tmp = Hap5.num{1}; 213 save('Spring5Coeff.txt','tmp','-ascii'); 214 break; 215end 216 217%8th order linkwitz riley x-over, used for chirp straightening 218[b a] = butter(4,Fc/2/Nyquist); 219Hxl = tf(conv(b,b),conv(a,a),Ts); 220[b a] = butter(4,Fc/2/Nyquist,'high'); 221Hxh = tf(conv(b,b),conv(a,a),Ts); 222 223%anti-aliasing filter 224%[b a] = ellip(10,1.5,80,Fc/Nyquist);%this needs to cutoff at Fc for each spring 225[b a] = cheby1(10,1.5,Fc/Nyquist); 226Haa = tf(b,a,Ts); 227%freqz(b,a,Fs); 228%break 229%we won't model the downsampling for this exercize 230 231if(0) 232 [y t] = impulse(Haa,.06); 233 %[y t] = lsim(Hxh,y,t); 234 [y t] = lsim(Hap,y,t); 235 specgram(y,256,Fs,hanning(256),256-32) 236 break; 237end 238 239%chirp straightening 240wxo = pi*(Fc/2)/Nyquist; %crossover frequency (in range [0,1]*pi) 241Da = k*M*(1-al^2)/(1+2*al*cos(wxo*k)+al^2);%delay of allpass filter at Xover freq 242Dad = Da - floor(Da);%decimal part of Da 243astrt = (-Dad*cos(wxo) + sqrt(1-Dad^2+Dad^2*cos(wxo)^2))/(1+Dad); %fractional delay line coeff. 244Hstrt = tf([astrt 1],[1 astrt zeros(1,floor(Da)-1)],Ts);%tf([astrt 1],[1 astrt zeros(1,k*M)],Ts); 245 246if(0) 247 [y t] = impulse(Haa,.06); 248 [yh t] = lsim(Hxh,y,t); 249 [yh t] = lsim(Hap,yh,t); 250 [yl t] = lsim(Hxl,y,t); 251 [yl t] = lsim(Hstrt,yl,t); 252 y = yl+yh; 253 specgram(y,256,Fs,hanning(256),256-32) 254 break; 255end 256 257%feedback delay 258Dlf = ceil(.056*Fs);%low frequency FB delay in samples (~56ms) 259%spring 4,5 will be 74,91ms 260glf = -.64; 261Hd = tf([glf],[1 zeros(1,Dlf)],Ts);%pure delay 262 263%putting it all together 264if(0) 265 Gap = Hxh*Hap + Hxl*Hstrt; 266 G = Gap/(1+Hd*Gap); 267 %bode(G) 268 %pause; 269 %[y t] = impulse(G,1,.6); 270 %[y t] = impulse(1/(1+Hd),.5); 271 %specgram(y,256,Fs,hanning(256),256-32) 272 %break; 273 [y t] = impulse(Haa,.5); 274 [y1 t] = lsim(Gap,y,t); 275 [y2 t] = lsim(1/(1+Hd*Gap),y1,t); 276 specgram(y2,256,Fs,hanning(256),256-32) 277 %[y t] = impulse(Gap,1,.5); 278 %specgram(y,[],44100); 279 %specgram(cf,512,44100,hanning(512),496) 280end 281 282%eq the LF chirp/anti-imaging filter 283[b a] = ellip(1,1.5,50,[180/Nyquist*4 320/Nyquist*4]); 284[b a] = ellip(1,1,50,[200/Nyquist*4 4000/Nyquist*4]); 285%900 3400 286Hai = tf(b,a,Fs/4); 287b2 = [b(1) 0 0 0 b(2) 0 0 0 b(3)]; 288a2 = [a(1) 0 0 0 a(2) 0 0 0 a(3)]; 289 290%freqz(b2,a2,Fs/4); 291%bode(Hai) 292%y = impulse(Hai); 293%specgram(y,64,Fs/4,hanning(64),48) 294 295%I think the AI filter should be another cheby1 296 297%the high frequency side 298%[b a] = ellip(5,3,6,[Nyquist/Fc Nyquist/upB]);a% 299al = -.64; 300k = 2; 301M = 200; 302Hap = (tf([al zeros(1,k-1) 1],[1 zeros(1,k-1) al],Ts))^M;%old way 303%[y t] = impulse(Hap,.05); 304%specgram(y,256,Fs,hanning(256),256-32); 305 306 307disp "Starting over" 308 309%new design decimate x4 310%5500hz k=2; al=-.7; M=19; D=56ms 311%5500hz k=2; al=-.5; M=24; D=88ms 312%5500hz k=2; al=-.3; M=21; D=63ms 313%3700hz k=3; al=-3.6;M=42; D=63.2ms 314%2800hz k=4; al=-.6; M=22; D=63.15ms 315nspring = 5; 316ks = [2 2 2 3 4]; 317als = [.7 .5 .3 3.6 .6]; 318 319%low side 320k=3; al=-.36; M=42; D=.063; 321z = roots([al zeros(1,k-1) 1]); 322p = roots([1 zeros(1,k-1) al]); 323b = poly(repmat(z',1,M)); 324a = poly(repmat(p',1,M)); 325Hap = zpk(repmat(z,M,1),repmat(p,M,1),1,Ts*4); 326%Hap = (tf(b,a,Ts*4));%*tf([1],[1 zeros(1,D*Fs/4)],Ts*4);%delay 327Hap2 = (tf([al zeros(1,k-1) 1],[1 zeros(1,k-1) al],Ts*4))^M;%old way 328[y t] = impulse(Hap,.06); 329specgram(y,64,Fs/4,hanning(64),64-8) 330 331 332k=2; al=-.7; M=19; D=.056; 333z = roots([al zeros(1,k-1) 1]); 334p = roots([1 zeros(1,k-1) al]); 335b = poly(repmat(z',1,M)); 336a = poly(repmat(p',1,M)); 337%Hap2 = (tf(b,a,Ts*4))*tf([1],[1 zeros(1,D*Fs/4)],Ts*4); 338 339figure(2) 340[y2 t] = impulse(Hap2,.10); 341specgram(y2,64,Fs/4,hanning(64),64-8) 342figure(1) 343%Hap.den{1} - Hap2.den{1} 344%max(Hap.den{1} - Hap2.den{1}) 345 346%high side 347 348