-
Notifications
You must be signed in to change notification settings - Fork 1
/
karplus_strong.m
91 lines (73 loc) · 1.94 KB
/
karplus_strong.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%karplus_strong.m
%Program author: Akihiro Inui
%Program details: Implement Karplus-Strong Algorithm
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 1. Preamble and variable declaration
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Sample rate
Fs = 44100;
% Sample period
Ts = 1/Fs;
% Fundamental frequency
f0 = 110;
% Length of delay line
N = round(Fs/f0 - 0.5);
% Loss factor
rho = 0.95;
% Dynamics filter coefficient
R = 0.95;
% duration of simulation(sec)
tEnd = 3;
% length of output vector
M = tEnd*Fs;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 2. Create input vector
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create white noise input
v = -1 + (1+1)*rand(M,1);
% Create filtered vector u
u = zeros(N,1);
% Initialise output
u(1) = (1-R)*v(1);
% Calculate dynamics filter
for n = 2:N
u(n) = (1-R)*v(n) + R*u(n-1);
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 3. Create input/output vector
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create input vector
x = zeros(N,1);
% Create output vector
y = zeros(N+M,1);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 4. Calculate Karplus
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The equation of Karplus-Strong is %
% y(n) = x(n)+(rho/2)*(y(n-N)+y(n-(N+1))) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implement Karplus-Strong Algorithm
% for n < N, x(n)=u(n) and for n >= N, x(n)= 0
% n < N
for n = 1:N-1
x(n) = u(n);
y(n) = x(n);
end
% n = N, N+1
y(N) = (rho/2)*y(1);
y(N+1) = y(N);
%n >= N
for n = N+2:N+M
y(n) = rho*(y(n-N)+y(n-(N+1)))/2;
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 5. Listen to the output signal
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Normalise output
MaxAmp = max(abs(y));
y = y/MaxAmp;
% Sound
soundsc(y,Fs)