% State space model of non-homogeneous Drum head - Tabla, By 'Yasser'Seekings.

% There are three basic hits represented here :- Tin, a hit to the middle, Edge

% hit to the edge(!), and dmp which is a damped hit, on the real thing this

% exploits the higher modes of resonance but not here, any ideas?

n=13; % number of 'masses' 9 is the minimum

% for reasonable results

T=1/3; % 1/3 of a second spent on each 'hit'

fs=11000; % sample freq.

% Basic hit - tin

%----------------------------------------------------------------------------------------------------------

tin(1:2:2*n)=0.001*sin((0:n-1)*pi/n);

tin(2:2:2*n)=sin((0:n-1)*pi/(n-1));

% hit to the edge

%----------------------------------------------------------------------------------------------------------

edge(1:2:2*n)=0.001*sin((4:n+3)*pi/n);

edge(4)=[1];

edge(2*n)=0;

% damped hit

%----------------------------------------------------------------------------------------------------------

dmp(1:2*n)=zeros(2*n,1);

dmp(14)=[10];

% air damping factor - empirically

% derived

%----------------------------------------------------------------------------------------------------------

r=0.017e-3*ones(1,n);

% Masses of the system - empirically

% derived

%----------------------------------------------------------------------------------------------------------

ro1=4.56; % Mass per unit area of the sythi

ro2=0.226; % Mass per unit area of the lao

dr=0.13/n % width of the puri/number of masses

drs=dr^2;

tm=ones(1,7);

tm(1)=pi*drs*ro; % Calculates the mass of the middle mass

for i=2:4;tm(i)=pi*ro1*drs*((2*i)-1)/2;end; % Calculates the masses of one half of the synthi

for i=5:7;tm(i)=pi*ro2*drs*((2*i)-1)/2; end; % Calculates the masses of one half of the lao

for i=1:6; tmass(i)=tm(8-i);end; % Fills the tmass matrix, with the values calculated

for i=7:13; tmass(i)=tm(i-6);end; % above.

% overall tuning factor empirically derived.

% Equivelent to Tk, see text.

%----------------------------------------------------------------------------------------------------------

tuning=21000;

% spring constants of springs

% equal to the tuning

%----------------------------------------------------------------------------------------------------------

k(1:14)=tuning*ones(14,1);

% state space matrices a,b,c,d

%----------------------------------------------------------------------------------------------------------

a=zeros(2*n);

b=zeros(size(a(:,1)));

for i=1:2:2*n; b(i)=10; end;

c=eye(size(a)); d=b;

for j=1:2:2*n % Fills A matrix for odd-numbered states.

a(j,j+1)=1;

end

% first mass

a(2,[1:3])=[-(k(1)+k(2))/tmass(1) -r(1)/tmass(1) k(2)/tmass(1)];

% last mass

a(2*n,[2*n-3:2*n])=[k(n)/tmass(n) 0 -(k(n)+k(n+1))/tmass(n) -r(n)/tmass(n)]; %last mass

% middle masses

for j=4:2:2*n-2

i=j/2;

a(j,[j-3:j+1])=[k(i)/tmass(i) 0 -(k(i)+k(i+1))/tmass(i) -r(i)/tmass(i) k(i+1)/tmass(i)];

end

% A matrix da for damped

% hit note different value for r(4)

%----------------------------------------------------------------------------------------------------------

r(5)=0.9;

da=zeros(2*n);

for j=1:2:2*n

da(j,j+1)=1;

end

% first mass

da(2,[1:3])=[-(k(1)+k(2))/tmass(1) -r(1)/tmass(1) k(2)/tmass(1)];

% last mass

da(2*n,[2*n-3:2*n])=[k(n)/tmass(n) 0 -(k(n)+k(n+1))/tmass(n) -r(n)/tmass(n)]; %last mass

%middle masses

for j=4:2:2*n-2

i=j/2;

da(j,[j-3:j+1])=[k(i)/tmass(i) 0 -(k(i)+k(i+1))/tmass(i) -r(i)/tmass(i) k(i+1)/tmass(i)];

end

% the business bit

% this simulates the drum head

% note the last values of one hit

% consitute to the next hit giving

% a more realistic situation

%----------------------------------------------------------------------------------------------------------

t=0:1/fs:T;

ts=fs*T

y=initial(a,b,c,d,tin,t);

y2=initial(a,b,c,d,edge+y(ts,:),t);

y3=initial(da,b,c,d,dmp+y2(ts,:),t);

% here we add the velocities

% of all the masses at a given time

% instant to give the output tsound

%----------------------------------------------------------------------------------------------------------

tsound(1:ts+1)=sum(y(:,2:2:2*n)');

tsound(ts+2:2*ts+2)=sum(y2(:,2:2:2*n)');

tsound(2*ts+3:3*ts+3)=sum(y3(:,2:2:2*n)');

pause;

saxis('auto');

sound(tsound,fs);

function [a,b,c,d,n,cn,cood]=drum(noincircum, noinradius,tuning);

% Yasser 96.

% Function which creates, state-space matrices a,b,c,d and

% the order n for the dahina, the left hand tabla, a

% classical Indian drum. Left hand arguements are the

% number of masses in the circumference, the number of

% masses in the radius+1 (not including the middle one,

% and finally the tuning constant try 9000, for this.

% Oh and the a matrix is sparse.

%connectivity matrix- cn

n=(noincircum*noinradius)+1; % n= total number of masses

dr=0.13/((2*noinradius)+1); % dr=width of puri/number of

% masses in diameter

%Mass matrix

m=ones(n,1); % Initializes mass matrix of the correct size

ro1=4.56; % Kg/m^2 of the synthi

ro2=0.226; % Kg/m^2 of the lao

m(1)=ro1*(dr^2)*pi; % Calulates the mass of the centre mass

cood=zeros(n,2); % Intialializes the coordinate matrix

w=2*pi/noincircum; % add coodinates, w is angular rotatation

% This calculates the coordinates of each of

for k=1:noinradius % the masses

for i=2:noincircum+1

cood(i+(noincircum*(k-1)),:)=[k*dr*sin((i-1)*w) k*dr*cos((i-1)*w)];

end;

end;

cn=sparse(zeros(n+1,n)); % This Initializes a sparse matrix for the

% connectivity information (cn)

cn(2:noincircum+1,1)=ones(noincircum,1); % fills in the cn matrix for centre (1) mass

% this fills the cn matrix

for k=0:noinradius-1;

pt1=2+k*noincircum;

pt2=k*noincircum+1+noincircum;

for i=pt1:pt2;

if i==pt1 cn(i+noincircum-1,i)=1;

else cn(i-1,i)=1; end; %left mass

if i==pt2 cn(i-noincircum+1,i)=1;

else cn(i+1,i)=1; end; %right mass

if k==0 cn(1,i)=1; % inside mass

else cn(i-noincircum ,i)=1; end;

if k==noinradius-1 cn(n+1,i)=1; % boundary

else cn(i+noincircum ,i)=1; end; % or an outside mass

m(i)=((2*k)+3)*(pi*dr^2)/noincircum; % This calculates the area of the puri that % each mass is approximating then it is % multiplied by the correct mass per unit

% area, depending on weather

if k<(noinradius-1)/2 m(i)=m(i)*ro1; % the mass represents the sythi or lao

else m(i)=m(i)*ro2; end;

end;%for

end; %for

%spring constant

%tuning=11000; % This sets the tuning/spring % constant

% damping % This is the ammount of % damping for the model

r=0.01;

a=sparse(zeros(n)); % Inializes a sparse matrix for A

nofcm=sum(cn); % Add up the columns of cn to give the number of % masses or boundaries one particular mass is % connected to.

for j=2:2:2*n

i=j/2;

a(j-1,j)=1; % puts in the odd states

a(j,j)=-r/m(i); % puts in velocity term

a(j,j-1)=-nofcm(i)*tuning/m(i); % puts in acceration -the number % of mass connected to mass i

end;

[l,u]=find(cn(1:n,1:n)); % Finds all the ones in the cn % matrix excluding those relating % to boundary conditions, there % indices are returned in l and u

for j=1:length(l)

a(2*l(j),2*u(j)-1)=tuning/m(l(j)); % puts in the rest of the % displacement terms

end;

b=zeros(size(a(:,1))); % Initializes B matrix

c=zeros(size(a)); % Initializes C matrix

for i=1:2*n c(i,i)=[1]; end; % and selects the outputs

d=b; % Initializes D matrix

[a,b,c,d,n]=drum(8,4,25000); % Calls the drum function to make the state % matrices and n the number of masses

fs=11000; % Sample frequency

T=0.3; % Length of each response to each hit

ts=fs*T; % The number of samples in each response

t=[0:1/fs:T]; % The time vector

% Below are the initial conditions for three different hits

% Basic hit - tin

%-----------------------------------

tin=zeros(1,2*n);

tin(16)=[40];

% hit to the edge

%-----------------------------------

edge=zeros(1,2*n);

edge(40)=[40];

% hit damped hit

%-----------------------------------

dmp=zeros(1,2*n);

dmp(52)=[60];

% This bit calculate the systems response to each hit, first without ‘finger factor’ (see text)

y=initial(full(a),b,c,d,tin,t);

y2=initial(full(a),b,c,d,edge+y(ts,:),t);

Then with the ‘finger factor’ !!!

a(8,8)=a(8,8)*10000; % ‘finger factor’

y3=initial(full(a),b,c,d,tin+y2(ts,:),t);

y4=initial(full(a),b,c,d,dmp+y3(ts,:),t);

% As you can see the initial conditions of the new hit is added to last state of the previous hit.

y=sum(y');

y2=sum(y2');

y3=sum(y3');

y4=sum(y4');

y4=y4(1:1100);

%This arranges them nicely

phrase=[y y3 y2 y4];

n=1;

n=n*4; % number of elements, aways a multiple of 4.

len=1/n;

kd=(0.5/len).*[1 -1;-1 1]; % element stiffness matrix

skd=zeros(n+1);

% This assembles the element stiffness

% matrices into a global stiffness matrix

for i=1:n

skd(i,i)=kd(1,1)+skd(i,i);

skd(i,i+1)=kd(1,2)+skd(i,i+1);

skd(i+1,i)=kd(2,1)+skd(i+1,i);

skd(i+1,i+1)=kd(2,2)+skd(i+1,i+1);

end;

% the next bit creates the mass matrix

ro1=1;

ro2=1;

mass=zeros(n+1);

for i=1:n+1

if(i>(n/4)+1)&(i<(3*n/4)+1) mass(i,i)=ro1;

else mass(i,i)=ro2;

end;

end;

% add boundary conditions

Km=mass(2:n,2:n);

Kd=skd(2:n,2:n);

% solve for eigenvalues and vectors

[v,l]=eig(Kd,Km);

[ev]=diag(l);

% plot results

figure(1);

subplot(221)

plot([0 ;v(:,1); 0],'c');

axis([1 n+1.1 -1 1]);

%title(sprintf('Eigenvector associated with lamda=%g',ev(1)));

subplot(222)

plot([0 ;v(:,2); 0],'c');

axis([1 n+1.1 -1 1]);

%title(sprintf('Eigenvector associated with lamda=%g',ev(2)));

subplot(223)

plot([0 ;v(:,3); 0],'c');

axis([1 n+1.1 -1 1]);

%title(sprintf('Eigenvector associated with lamda=%g',ev(3)));