% E04_8          26apr17    02apr19  20apr21
% comb filters (FIR & IIR)
clear all;
return;

%%  basic comb filters - FIR
for i = 0:10;    % ex.: 0:10   note: i=L-1
    h = [1  zeros(1,i)  1]/2;  
% h = [1 zeros(1,i) 1 zeros(1,i) 1 zeros(1,i) 1]/4;
    Hm = abs(fft(h,512));
    subplot(221); zplane(h,1);
    subplot(222); plot(Hm(1:257)); axis('tight');
    h = [1  zeros(1,i)  -1]/2;  
% h = [1 zeros(1,i) -1 zeros(1,i) 1 zeros(1,i) -1]/4;
    Hm = abs(fft(h,512));
    subplot(223); zplane(h,1);
    subplot(224); plot(Hm(1:257)); axis('tight');
    pause();
end;

%%  basic comb filters - IIR
for i = 0:10;    % ex.: 0:10
    b = [1  zeros(1,i)  1];   a = [1  zeros(1,i)  -.7];  
    Hm = abs(freqz(b,a,256));
    subplot(221); zplane(b,a);
    subplot(222); plot(Hm); axis('tight');

    b = [1  zeros(1,i)  -1];   a = [1  zeros(1,i)  -.7];  
    Hm = abs(freqz(b,a,256));
    subplot(223); zplane(b,a);
    subplot(224); plot(Hm); axis('tight');
    pause();
end;

%%  echo - reverberation
[x,fs] = audioread('dsp01.wav');   % Fs = 44 kHz
x = x(1:4:end); fs = fs/4; % reduce computational load
zz = zeros(1,1000);
 h = [1 zz -1]; % test also last coeff. <1 or <0
% h = [1 zz 1 zz .7 zz zz .3];
y = filter(h,1,x);
soundsc(y,fs);

%%  filtering example
x = sin(.123*pi*[0:1:199]);
thr = .01;   x(x>=thr) = 1;  x(x<-thr) = -1;
x = x + randn(1,200)*.1; % try  .1 ... 1
X = abs(fft(x));
 alpha = .72;    beta = cos(.123*pi);  
 b = [1  0  -1] * (1-alpha) /2;   a = [1  -beta*(1+alpha)  alpha];   % IIR bandpass, 4(5)30
zz = zeros(1,7);
h1 = [1 zz -1]/2;  % basic highpass, L=7 --> peaks at pi/8 and multiples
h2 = [1 zz -1 zz 1 zz -1 zz 1 zz -1 zz 1 zz -1]/8;  % running avg, z to -z, L=7
H1 = abs(fft(h1,128));
H2 = abs(fft(h2,128));
figure(1);
subplot(411); plot(X(1:100)); axis('tight');
subplot(412); plot(abs(freqz(b,a)));  axis('tight');
subplot(413); plot(H1(1:64)); axis('tight');
subplot(414); plot(H2(1:64)); axis('tight');
y0 = filter(b,a,x);
y1 = filter(h1,1,x);
y2 = filter(h2,1,x);
figure(2);
subplot(411); plot(x);
subplot(412); plot(y0);
subplot(413); plot(y1);
subplot(414); plot(y2);


%%  IIR comb
%  Karplus-Strong "pizzicato" (plucked string) generator
clf;

Fs = 8192;  % note: fundamental freq. is close to Fs/i (lowest freq. pole at pi/(i/2) )
% for i = 20:5:40;  % ex.: 0:20
% for i=[19 17 15 14 12 11 10 9]; % poor musical scale: too low resolution (Fs can be used for fine tuning)
 for i = round(Fs ./ (261.6 / 4 * 2.^([0 2 4 5 7 9 11 12]/12)));   % scale starts from C 1 octave below 440 Hz A:  220*2^(3/12))
    x = [randn(1,i+2), zeros(1,9900)];
 %   x = [ones(1,1), zeros(1,9900)];
 %   x = [linspace(0,1,i), zeros(1,9900)];
    b = 1;  a = [1 zeros(1,i) -.5 -.5];
 %   h = .5;  b = 1;  a = [1 zeros(1,i) h-1 -h];
 %   k = .5;  h = .5;  b = [1 zeros(1,floor(i/2)) k-1 -k];  a = [1 zeros(1,i) h-1 -h];
    figure(1); clf; zplane(b,a); pause(.1);
    figure(2); freqz(b,a); pause(.1);
    figure(3); impz(b,a,500); pause(.1);
    y = filter(b,a,x);
    figure(4); plot(y);
    sound(y,Fs); pause(.5);  % note Fs default = 8192 Hz
    end
return;

% (TI_appl_note.pdf)
% Variant: excite the line with a noise pulse with variable rise
% and decay times. It makes some cool sounds, ranging from 
% quasi-bowed attacks, to a buzzy "jews harp" effects.
% It's also possible to generate sustained sounds, by continuing 
% to feed some noise into the system, and to lowpass filter the 
% noise, which alters the color of the sound.
%
% see also code available on Mathworks.com
