%% Audio effects
% Table of Contents
% Echo FIR & IIR (comb)
% Allpass
% Schroeder's reverb
% Flanger
% Chorus
% Tremolo
% WHA WHA filter
% Amplitude modulation
% Frequency modulation
%
%
%% Echo FIR & IIR (comb)
%
% Echo filters based on the delay line struct

%% FIR echo filter 
disp('FIR Echo filter of order (D) and gain (g):')
D = 10   % order
g = 0.4  % gain

% coefficients
a = 1;                    % AR coeff
b  = [1 zeros(1,D-1) g];  % MA coeff

% analysis
freqz(a,b,256)  % freq resp.
zplane(a,b)     % zeros/poles
impz(b,a,20)    % impulse resp.

%

%% IIR echo filter 
disp('IIR Echo filter of order (D) and gain (g):')
D = 10   % order
g = 0.8  % gain

% coefficients
a  = [1 zeros(1,D-1) g];  % AR coeff
b = 1;                    % MA coeff

% analysis
freqz(a,b,256)  % freq resp.
zplane(a,b)     % % zero-pole pattern 
impz(b,a,100)   % impulse resp.

%
%% Allpass
%

% allpass filter of order D

disp('Allpass filter of order (D) and gain (g):')
D = 6   % order
g = 0.8  % gain

% Set coeffs depending on g and D
b = [g zeros(1,D-1) 1];
a = [1 zeros(1,D-1) g];
        
impz(b,a,80)
zplane(b,a)    % zero-pole pattern

%filter the input signal
%y=filter(b,a,x);

%
%% Schroeder's reverb
% A reverberator based on Schroeder's design which consists of 4 parallel feedback comb filters in series with 2 allpass filters
%

%This is a reverberator based on Schroeder's design which consists of 4 parallel feedback comb filters 
%in series with 2 cascaded all pass filters.
%
%The structure is:  [y,b,a] = schroeder2(x,cg,cd,ag,ad,k)
%
%where x = the input signal
%      cg = a vector of length 4 which contains the gain of each of the comb filters (should be less than 1 for stability)	
%      cd = a vector of length 4 which contains the delay of each of the comb filters 
%      ag = the gain of the allpass filters (this should be less than 1 for stability)
%      ad = a vector of length 2 which contains the delay of each of the allpass filters 
%      k = the gain factor of the direct signal
%      y = the output signal
%      b = the numerator coefficients of the transfer function
%      a = the denominator coefficients of the transfer function
%
%
% Gautham J. Mysore - gauthamjm@yahoo.com
%

cg1 = [.75 .75 .75 .75 .75 .75];
cg = cg1; 
cd = [29 37 44 50 68 79];
ag = [.75 .75];
ad = [27 31];
k = 1;
x = [1 zeros(1,100)];

% send the input to each of the 4 comb filters separately
[outcomb1,b1,a1] = fbcomb(x,cg(1),cd(1));
[outcomb2,b2,a2] = fbcomb(x,cg(2),cd(2));
[outcomb3,b3,a3] = fbcomb(x,cg(3),cd(3));
[outcomb4,b4,a4] = fbcomb(x,cg(4),cd(4));

% sum the ouptut of the 4 comb filters
apinput = outcomb1 + outcomb2 + outcomb3 + outcomb4; 

%find the combined filter coefficients of the the comb filters
[b,a]=parallelcoefficients(b1,a1,b2,a2);
[b,a]=parallelcoefficients(b,a,b3,a3);
[b,a]=parallelcoefficients(b,a,b4,a4);

% send the output of the comb filters to the allpass filters
[y,b5,a5] = allpass(apinput,ag,ad(1));
[y,b6,a6] = allpass(y,ag,ad(2));

%find the combined filter coefficients of the the comb filters in series with the allpass filters
[b,a]=seriescoefficients(b,a,b5,a5);
[b,a]=seriescoefficients(b,a,b6,a6);

% add the scaled direct signal
y = y + k*x;

% normalize the output signal
y = y/max(y);

%
%% Flanger
%
% Creates a single delay with the delay time oscillating from either 0-10 ms or 0-15 ms at 0.1 - 5 Hz this is not synthesisable unless buffering is used
% do calculations with sampling frequency to convert delay in samples into miliseconds (need some 44.1kHz samples)
%

% load file
filename = 'music/acoustic.wav'
[x,fs] = audioread(filename);
info = audioinfo(filename)
plot(x)
xlabel('Sample Number')
ylabel('Amplitude')
title('acoustic.wav')

% current sample is 11kHz so 0-10 ms is 0 - 110 samples
max_time_delay = 0.01;   % 10 ms max delay in seconds
Fd = 1;                  % rate of flange in Hz
N = length(x);           % tot samples

% convert delay in ms to max delay in samples
D = round(max_time_delay*fs); 

% low frequency oscillating delay (cast to int)
n = 1:N; % indexes
d = floor(D/2*(1-cos(2*pi*n*(Fd/fs))));

% create empty out vector
y = zeros(1,N);      

% to avoid referencing of negative samples
y(1:D) = x(1:D); 

amp = 0.7; % suggested coefficient from page 71 DAFX

% for each sample
for i = (D+1):N
    % add delayed sample
    y(i) =  x(i) + amp * x(i-d(i));   
end

% plots
tms = (1:N)/fs; 
plot(tms,x,'r',tms,y,'b');
title('Flanger and original Signal');

%

% original
sound(x,fs)

%

% flanger
sound(y,fs)

%
%% Chorus
%
% Chorus is an effect that simulates the presence of several sources playing in imperfect unison. It is implemented by combining the original signal with several copies delayed by a randomly varied number of samples. Typical delays are in the range D = 10 - 25 ms, and they are modified by small, slowly varying random variations.
% ﻿For the linear interpolation case, we need to keep track of two successive random values R, say r\left\lbrack i\right\rbrack and r\left\lbrack i+1\right\rbrack, and connect them linearly
% Because the slope of the straight line between r\left\lbrack i\right\rbrack and r\left\lbrack i+1\right\rbrack is \left(r\left\lbrack i\right\rbrack +r\left\lbrack i+1\right\rbrack \right)/D, the linearly interpolated samples will be
%     d\left\lbrack i\right\rbrack =r\left\lbrack i\right\rbrack +\left(r\left\lbrack i+1\right\rbrack -r\left\lbrack i\right\rbrack \right)\frac{t}{D},\;\;\;t=0,1,\ldotp \ldotp \ldotp ,D-1
%

% load file
filename = 'music/acoustic.wav'
[x,fs] = audioread(filename);
info = audioinfo(filename)
plot(x)
xlabel('Sample Number')
ylabel('Amplitude')
title('acoustic.wav')

% current sample is 11kHz
F = 0.02;          % rate of variation in [cycl/samp]
D = floor(1/F);     % period of random generator
N = length(x)- mod(length(x),D);  % tot samples (multip. of D)
x = x(1:N);   % tot samples is multiple of D

% random variations for 2 sounds
R = N * F;
r1 = rand(1,R)-0.5; % zero-mean random numbers
r2 = rand(1,R)-0.5; % zero-mean random numbers

% create 2 empty delays vectors
d1 = zeros(1,N); 
d2 = zeros(1,N);

% random delays using linear interpolation in between
for i = 1:R-1
    % sound first delays
    line = r1(i)+(r1(i+1)-r1(i))*(0:D-1)/D;
    d1((i-1)*D+1:i*D) = round(D*(0.5 + line));
    
    % sound second delays
    line = r2(i)+(r2(i+1)-r2(i))*(0:D-1)/D;
    d2((i-1)*D+1:i*D) = round(D*(0.5 + line));
end

% create empty out vector
y = zeros(N,1);      

% to avoid referencing of negative samples
y(1:D) = x(1:D); 

amp = 0.7; % suggested coefficient from page 71 DAFX

% for each sample
for i = (D+1):N
    % first random delayed sound
    y1 = x(i-d1(i));
    
    % second random delayed sound
    y2 = x(i-d2(i));
    
    % sum up
    y(i) = x(i) + amp/4*(y1+y2);   % add delayed sample
end

% plots
tms = (1:N)/fs; 
plot(tms,x,'r',tms,y,'b');
title('Flanger and original Signal');

%

% original
sound(x,fs)

%

% chorus
sound(y,fs)

%
%% Tremolo
%

Fc = 5;
alpha = 0.5;
trem=(1+ alpha*sin(2*pi*n*(Fc/fs)))';
y = trem.*x;
% write output
audiowrite(outfile,y,fs);
figure(1)
plot(x,'r');
title('Original Signal');
figure(2)
hold on
plot(x,'r');
plot(y,'b');
title('Tremolo and original Signal');

%

soundsc(x) % sound tremolo

%
%% WHA WHA filter
% BP filter with narrow pass band, Fc oscillates up and down the spectrum. Difference equation taken from DAFX chapter 2.
%

% wah_wah.m   state variable band pass
% Changing this from a BP to a BS (notch instead of a bandpass) converts 
% this effect to a phaser
%
% yl(n) = F1*yb(n) + yl(n-1)
% yb(n) = F1*yh(n) + yb(n-1)
% yh(n) = x(n) - yl(n-1) - Q1*yb(n-1)
%
% vary Fc from 500 to 5000 Hz
        % 44100 samples per sec

infile = 'acoustic.wav';

% read in wav sample
[x,Fs] = audioread(infile);

%%%%%%% EFFECT COEFFICIENTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lower the damping factor the smaller the pass band
damp = 0.05;

% min and max centre cutoff frequency of variable bandpass filter
minf=500;
maxf=3000;

% wah frequency, how many Hz per second are cycled through
Fw = 2000; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% change in centre frequency per sample (Hz)
delta = Fw/Fs;
%0.1 => at 44100 samples per second should mean  4.41kHz Fc shift per sec

% create triangle wave of centre frequency values
Fc = minf:delta:maxf;
while length(Fc) < length(x)
    Fc = [ Fc (maxf:-delta:minf) ];
    Fc = [ Fc (minf:delta:maxf) ];
end

% trim tri wave to size of input
Fc = Fc(1:length(x));

% difference equation coefficients
F1 = 2*sin((pi*Fc(1))/Fs);  % must be recalculated each time Fc changes
Q1 = 2*damp;                % this dictates size of the pass bands

yh=zeros(size(x));          % create emptly out vectors
yb=zeros(size(x));
yl=zeros(size(x));

% first sample, to avoid referencing of negative signals
yh(1) = x(1);
yb(1) = F1*yh(1);
yl(1) = F1*yb(1);

% apply difference equation to the sample
for n=2:length(x)
    yh(n) = x(n) - yl(n-1) - Q1*yb(n-1);
    yb(n) = F1*yh(n) + yb(n-1);
    yl(n) = F1*yb(n) + yl(n-1);
    F1 = 2*sin((pi*Fc(n))/Fs);
end

%normalise
maxyb = max(abs(yb));
yb = yb/maxyb;

% write output wav files
%audiowrite('out_wah.wav',yb,Fs);

figure(1)
hold on
plot(x);
plot(yb);
title('Wah-wah and original Signal');

%
%% Amplitude modulation
%

N = 2000;
Om_e = 2/N; % Omega env freq
Om = 0.01;  % Omega fundamental freq
n = 0:N-1;
A = 1 + .25*cos(2*pi*Om_e*n); % oscillator
x = sin(2*pi*Om*n);           % signal
y = A.*x;
plot(y)
hold on
plot(A,'-r')
plot(-A,'-r')
hold off
figure(2)
clf(2)
Om_e = 3.5/N;
A = 1 + .25*cos(2*pi*Om_e*n); % oscillator
y = A.*x;
plot(y)
hold on
plot(A,'-r')
plot(-A,'-r')
hold off

%
%% Frequency modulation
%

N = 1000;
Om_c = 0.02;    % carrier freq
Am = 0.5*Om_c;  % Amplitude
Om_m = 0.004;   % frequency
n = 0:N-1;
Om = Om_c + Am*sin(2*pi*Om_m*n); % oscillator
y = sin(2*pi*Om.*n);           % signal
plot(y)
hold on
plot(y,'-b')
plot(Om/Om_c+0.5,'--r','LineWidth',2)

figure(1)
clf(1)
N = 1000;
Om_c = 0.02;    % carrier freq
Am = 0.5*Om_c;  % Amplitude
Om_m = 0.004;   % frequency
n = 0:N-1;
Om = Om_c + Am*ceil(sin(2*pi*Om_m*n)); % oscillator
y = sin(2*pi*Om.*n);           % signal
plot(y)
hold on
plot(y,'-b')
plot(Om/Om_c,'--r','LineWidth',2)