Hi everyone,
I am trying to make my own simulating tool using GNU Octave and transfer matrices methode. I am encountering difficulties in validating my simple model (infinite baffle) because it does not match hornresp results. i think i am not that far from a correct model but not here yet. as you can see phase and SPL are not the same as hornresp.
below is the code with some comments. Note that i do not think problem come from my radiatoin impedance since i used this function for analytical modeling and results mached hornresp with my analytical model.
I am trying to make my own simulating tool using GNU Octave and transfer matrices methode. I am encountering difficulties in validating my simple model (infinite baffle) because it does not match hornresp results. i think i am not that far from a correct model but not here yet. as you can see phase and SPL are not the same as hornresp.
below is the code with some comments. Note that i do not think problem come from my radiatoin impedance since i used this function for analytical modeling and results mached hornresp with my analytical model.
Code:
clear all; close all; clc;
[c, rho, P, P_ref, Air_Imp, f, w] = Generals(); % initialise constantes and frequenies
wc = w / c; % wave number
driver = 'RCF_MB15N405.csv'; % Selecting driver
r = 1; % Distance to from the source
eg = 2.83; % Input voltage
% Thiele/Small parameters
[Sd, fs, Vas, Qms, Bl, Re, Le, Xmax, Mms, Qts, Qes, Rms, n, dr] = TS_extract(driver, 1);
Cms = Vas / (rho * c^2 * Sd^2) % Calculating mechanical compliance
Pg=eg*Bl/(Sd*Re);
mas=rho*Sd;
% Radiation impedance uing bessel and struve functions choose
[zar,zmr]=RadImpC(Sd);
% Loop for each frequency
for k = 1:length(w)
% Transfer matrices for sub domains
tre(:, :, k) = [1, Re; 0, 1]; % voice coil resistance
tle(:, :, k) = [1, 1i*w(k)*Le; 0, 1]; % Inductance
tem(:, :, k) = [0, Bl; 1/Bl, 0]; % gyrator
tmm(:, :, k) = [1, 1i*w(k)*Mms; 0, 1]; % Mobil mass
tcm(:, :, k) = [1, 1/(1i*w(k)*Cms); 0, 1]; % Compliance
trm(:, :, k) = [1, Rms; 0, 1]; % mechanic losses
% Electromechanic chain
tdr(:, :, k) = tre(:, :, k) * tem(:, :, k) * tmm(:, :, k) * tcm(:, :, k) * trm(:, :, k);
% Calculating Iin [A] electrical current from the source from electro-mechanic
% matrix [A,B ; C,D]
t12 = tdr(1,2,k); % A
t22 = tdr(2,2,k); % B
zin = t12 / t22; % Impedance seen from source
ein = eg; % Input voltage
iin = ein / zin; % Input current
In = [ein; iin];
% Acoustic domain
tma(:, :, k) = [Sd, 0; 0, 1/Sd]; % Tranformer
tza(:, :, k) = [1, zar(k); 0, 1]; % Radiation Impedance
% Complete chain from eletrical domain to acoustical domain
Tt(:, :, k) = tre(:, :, k) * tem(:, :, k) * tmm(:, :, k) * tcm(:, :, k) * trm(:, :, k) * tma(:, :, k) * tza(:, :, k);
Ttinv = inv(Tt(:, :, k)); % Inverting matrix
Pu = Ttinv * In; % Calculating results
% Pression (pout) et vitesse acoustique (uout)
pout(k) = Pu(1); % Acoustic pressure at Sd
uout(k) = Pu(2); % Acoustic volume velocity
end
Pa = pout ./ (4 * pi * r); % Acoustic pressure at 1 [m] from the Sd
SPL = 20 * log10(abs(Pa)/ 20e-6); % Decibel SPL
G=1i.*w.*mas.*uout./Pg; % system gain
Phase_deg = unwrap(angle(G)) * 180/pi; % Phase
[fh,o,ph]=hornrespcsv("ibn4050.csv");
figure(1)
semilogx(f, SPL, 'LineWidth', 1.5, fh,o,'LineWidth',1.5);
xlabel('Frequencies (Hz)');
ylabel('SPL');
title(['SPL response']);
grid on;
figure(2)
semilogx(f, Phase_deg, 'LineWidth', 1.5,fh,ph,'LineWidth', 1.5);
xlabel('Frequencies (Hz)');
ylabel('Phase');
title(['Phase']);
grid on;
Last edited:
You are asking for someone to debug an octave script? I have never used octave and I suspect neither have a fair few others.
It's simpler to work with impedances which is in pretty much every acoustics text book. Might be quicker to implement that and compare against your octave script? Possibly more reliable since you will know exactly what assumptions are being made.
I don't know what is in hornresp and I think it is released as a windows binary(?) so you would have to ask the author who is active on the site.
Not much direct help I am afraid but hopefully it might help a bit with why you are not getting responses.
It's simpler to work with impedances which is in pretty much every acoustics text book. Might be quicker to implement that and compare against your octave script? Possibly more reliable since you will know exactly what assumptions are being made.
I don't know what is in hornresp and I think it is released as a windows binary(?) so you would have to ask the author who is active on the site.
Not much direct help I am afraid but hopefully it might help a bit with why you are not getting responses.
Parallelpen, did you try if Large Language Model ( AI ) could help?
Out of curiosity tried whether Grok understands any of it, and it spit out following.
I cannot comprehend what you do there so no idea whether this worked or not, but the result seems plausible so thought to post.
Problems in the code
After analyzing your code, I noticed the following key issues:
Out of curiosity tried whether Grok understands any of it, and it spit out following.
I cannot comprehend what you do there so no idea whether this worked or not, but the result seems plausible so thought to post.
Problems in the code
After analyzing your code, I noticed the following key issues:
- Missing voice coil inductance in the transfer matrix chain
You have defined the voice coil inductance matrix tle = [1, 1i*w(k)*Le; 0, 1], but you do not use it in any chain (tdr or Tt). The voice coil inductance (Le) is an important part of the speaker's electrical impedance, and omitting it affects the results, especially at higher frequencies. - Incorrect calculation of input impedance and current
You calculated the input impedance zin = t12 / t22, where t12 = tdr(1,2,k) and t22 = tdr(2,2,k), using only a partial chain tdr (the electro-mechanical part without the acoustic load). This does not take into account the entire system's load, including the radiation impedance (zar). Furthermore, the input current iin = ein / zin is based on this incorrect impedance, which distorts subsequent calculations. - Incorrect calculation of far-field pressure
You calculated the acoustic pressure at distance r using the formula Pa = pout ./ (4 * pi * r), where pout is the pressure on the membrane surface (pout = zar * uout). However, this does not correspond to the far-field pressure in an infinite baffle. The correct far-field pressure should be calculated from the volume velocity (uout) and not from the membrane pressure, and there is an error in the coefficient (you use 4 * pi * r, when it should be 2 * pi * r for the infinite baffle case). - Uncertainty in phase calculation
The phase is calculated using the formula G = 1i.*w.*mas.*uout./Pg, where mas = rho*Sd and Pg = eg*Bl/(Sd*Re). This seems arbitrary and may not correspond to Hornresp's calculation method. You probably want to calculate the phase of the far-field pressure directly.
Code:
clear all; close all; clc;
[c, rho, P, P_ref, Air_Imp, f, w] = Generals(); % Initialise constants and frequencies
wc = w / c; % Wave number
driver = 'RCF_MB15N405.csv'; % Selecting driver
r = 1; % Distance from the source
eg = 2.83; % Input voltage
% Thiele/Small parameters
[Sd, fs, Vas, Qms, Bl, Re, Le, Xmax, Mms, Qts, Qes, Rms, n, dr] = TS_extract(driver, 1);
Cms = Vas / (rho * c^2 * Sd^2); % Calculating mechanical compliance
% Radiation impedance using Bessel and Struve functions
[zar, zmr] = RadImpC(Sd);
% Loop for each frequency
for k = 1:length(w)
% Transfer matrices for subdomains
tre(:, :, k) = [1, Re; 0, 1]; % Voice coil resistance
tle(:, :, k) = [1, 1i*w(k)*Le; 0, 1]; % Inductance
tem(:, :, k) = [0, Bl; 1/Bl, 0]; % Gyrator
tmm(:, :, k) = [1, 1i*w(k)*Mms; 0, 1]; % Moving mass
tcm(:, :, k) = [1, 1/(1i*w(k)*Cms); 0, 1]; % Compliance
trm(:, :, k) = [1, Rms; 0, 1]; % Mechanical losses
tma(:, :, k) = [Sd, 0; 0, 1/Sd]; % Transformer
tza(:, :, k) = [1, zar(k); 0, 1]; % Radiation impedance
% Complete chain including inductance
Tt(:, :, k) = tle(:, :, k) * tre(:, :, k) * tem(:, :, k) * tmm(:, :, k) * tcm(:, :, k) * trm(:, :, k) * tma(:, :, k) * tza(:, :, k);
% Calculate volume velocity
u_a(k) = eg / (Tt(1,1,k) * zar(k) + Tt(1,2,k));
% Far-field pressure at distance r
p_far(k) = (1i * w(k) * rho / (2 * pi * r)) * u_a(k);
Pa(k) = (w(k) * rho / (2 * pi * r)) * abs(u_a(k)); % Magnitude for SPL
end
% Calculate SPL and phase
SPL = 20 * log10(abs(Pa) / 20e-6);
Phase_deg = unwrap(angle(p_far)) * 180 / pi;
% Load Hornresp data for comparison
[fh, o, ph] = hornrespcsv("ibn4050.csv");
% Plot SPL
figure(1)
semilogx(f, SPL, 'LineWidth', 1.5, fh, o, 'LineWidth', 1.5);
xlabel('Frequencies (Hz)');
ylabel('SPL (dB)');
title('SPL Response');
grid on;
legend('Octave Model', 'Hornresp');
% Plot phase
figure(2)
semilogx(f, Phase_deg, 'LineWidth', 1.5, fh, ph, 'LineWidth', 1.5);
xlabel('Frequencies (Hz)');
ylabel('Phase (degrees)');
title('Phase Response');
grid on;
legend('Octave Model', 'Hornresp');thank you for answering,
I was not clear in my first post. i am not asking for someone to debug my script but to explain to me how it should be done, like what is the correct way to calculate the input current, matrices order of multiplication and things like that.
octave is the same as matlab so i was hoping some people could be able to understand my code wich is here to show how i understand the concepts, i am not waiting for someone to send me a corrected script but someone to explain to me what is the correct way to use the matrix method for loudspeaker systems.
I wil try to explain what i am doing without codes maybe it is better.
I am trying to implement a simple model loudspeaker in an infinite baffle using trassmission matrices. I have to implement it correctly so that it represent accurately the physics of a loudspeaker in all three domains (electrical, mechanical and acoustical).
In a more concrete expression i want to use the following relationship to compute output pressure and volume velocity :
[pressure ; volume velocity] = Tinverted * [voltage ; curent]
So i need to know the correct voltage (here it is just a constant input voltage eg), the input current wich i need to calculate accurately and the inverted matrix reprensenting my complete model Tinverted.
. Electrical domain :
There for i am computing the electrical matrix buy multiplying the restance with the inductance matrix.
Telec = TRe * Tle
with TRe and Tle the resistance matrix and the inductance matrix created like stated in bjorn kolbrek book for impedance matrices so Tz=[1,Z ; 0,1]
In the code i sent i am not adding the inductance Le but this is not where i am mistaken, adding Le as stated above with does not make my results correct.
. Mechanical domain :
I am using the resistance (Z=Rms) giving me the mechanical resistance matrix TRms, the mass impedance (Z=j*w*Mms) giving me the matrix for the mechanical mass TMms and the compliance (Z=-j/w*Cms) giving me the matrix for the mechanical compliance TCms.
then i multiply those matrices in the following order : Tmech = TMms * TCms * TRms
My complete electromechanical driver matrix is then : Tdriver=Telec * Tem * Tmech
Where Tem is the electricomechanical coupling matrix Tem=[0,Bl ; 1/Bl,0] also taken from Kolbrek book.
. Calculating input current I :
Once this is done i am assuming i have a correct model of the electrical and mechanical sides so i calculate the input current. To compute the input current i am using the transmission matrix propriety wich is giving me the electrical impedance shorted output by using the driver matrix Tdriver components :
Zin = t12 / t22 where t12 and t22 are the Tdriver matrix components (Tdriver can be writen as : Tdriver = [t11,t12 ; t21,t22])
Once i have Zin i can calculate the input current :
I = eg / Zin where eg is the input voltage and Zin the previously calculated impedance seen from the source.
. Acoustical domain :
Now the last component missing is the inverted matrix Tinverted wich is just the inverse matrix of my complete system in all 3 domains. So i need to add the acoustical domain to my Tdriver matrix to create the complete system's matrix T :
T = Tdriver * Tma * Ta
where Ta is the matrix for the acoustical side, i am using the same method as for electrical and mechanical domain expect here Z is a radiation impedance calculated from bjorn kolbrek book section 14.4 for those who are curious.
Tma is the mechanoacoustical coupling matrix Tma = [Sd,0 ; 0,1/Sd]
. Final Computation :
so now i have a correct input voltage : input voltage = eg for all frequencies, a calculated input current I, and a inverted matrix Tinverted = invert of T.
[pressure ; volume velocity] = Tinverted * [voltage ; curent]
My problem is by doing so i do not have the correct output pressure nor the correct output volume velocity comparing to my analytical model nor the correct phase and spl compared to hornresp.
I hope those precisions can give you a more clear picture of my problem and my current understandings.
Last edited:
Thank you for answering,
1. The first point is not one of importance here i was voluntarly neglecting Le in the code i sent.
2. That is one thing i am not sure about, but what i can say is when calculating impedance Zin with the complete chain matrix T i don't get anything close from reality but noise. maybe i am not doing it the right way, an expert explanation would be really appreciated here.
3. In fact when using volume velocity (uout in my code) to compute the pressure my curve matches hornresp for the low frequencies but not for the active frequency range (see figure below). i decided to use pout in my code because i am assuming this is the way but again an expert explanation here about why the output pressure i am computing is not the correct one and how i should do it would be great.
4. I also tried but since my acoustic pressure is not correct it is not a solution for me, maybe when i will compute correct pressure and volume velocity. here also an expert poit of view can help me.
Last edited:
What do you mean by analytical model? Have you written down something like:
Zr = rfac*radimp(2*ka)Zc = complex(Rc, w*mc - sc/w)ZE = complex(RE, w*LE)Zm = Zr + ZcZM = phi2/ZmZI = ZE + ZMfrom an equivalent circuit, looking in an acoustics text book, or whatever. Put it in your frequency loop and evaluated the terms you wish to compare with the equivalent ones in your model?
my analytical model is another Octave code i did from reading Small paper about direct radiator enclosure and a few other books, using equivalent circuit and classical solving as if i was in school. This analytical model matches hornresp 100% phase wise and spl wise so i use this model for in depth comparisons.
Is there a problem printing out intermediate values in your two octave scripts in order to identify where the discrepency lies?
thanks for advice, I think i have solved my problem simply by calculating inpput current analyticaly as follow :
Total impedance =ze + ((Bl^2) / (zm + za * (Sd^2)))
input current = eg / zt;
But it would be cleaner if i can be able to compute the input current from my matrices. Any advice to do it properly ? I dont know if it is a normal result but the acoustic pressure i get as result from multiplying my matrices is just noise but the volume velocity is correct so i compute all acoustic results from volume velocity and all plots (SPL, phase and volume velocity compared to my analytical model that matches hornresp) are matching hornresp. If someone in the known can explain me what is going on here it would be very appreciated.
Total impedance =ze + ((Bl^2) / (zm + za * (Sd^2)))
input current = eg / zt;
But it would be cleaner if i can be able to compute the input current from my matrices. Any advice to do it properly ? I dont know if it is a normal result but the acoustic pressure i get as result from multiplying my matrices is just noise but the volume velocity is correct so i compute all acoustic results from volume velocity and all plots (SPL, phase and volume velocity compared to my analytical model that matches hornresp) are matching hornresp. If someone in the known can explain me what is going on here it would be very appreciated.
Last edited:
In the transfer matrix statement:
tmm( :, :, k) = [1, 1i*w(k)*Mms; 0, 1];
Mms should be Mmd
Mmd = Driver diaphragm and voice coil dynamic mechanical mass
Mmd equals Mms (total moving mass) minus air load
tmm( :, :, k) = [1, 1i*w(k)*Mms; 0, 1];
Mms should be Mmd
Mmd = Driver diaphragm and voice coil dynamic mechanical mass
Mmd equals Mms (total moving mass) minus air load
Thank you for answering, i rectified it in my code.
I think i finnally did the job. I don't know if my method is 100% accurate but results are 100% matching hornresp if calculating outputs from volume velocity uout and not using acoustical pressure pout.
What i am doing :
my matrix chain is the same as before, what changed is i am not calculating input current from the electromechanical matrix and the Zin=t12/t22 method but calculating directly by forcing an arbitrary output volume velocity and pressure and using [input voltage ; input current] = T * [pressure ; volume velocity] to compute virtuals inputs then calculating impedance seen from source Zin = input voltage / input current and then just calculating the real input current as
iin = eg / Zin.
With this input current iin i can compute a correct volume velocity uout without needing any analytical expression of acoustic load impedance witch i needed before when using analytical expression for Zin Total impedance =ze + ((Bl^2) / (zm + za * (Sd^2))).
Every thing is working great since i can compute everything from volume velocity but i still struggle to understand why acoustic pressure pout from matrix computation is incorrect but i guess it is a normal result since i am cheating a bit to compute input current iin, this must be the draw back of my method witch is not a big deal to be fair. Anyway if someone here knows how to compute input current without any analytical expressions and giving a correct pout and uout, let me know.
Now i can compute SPL, phase, excursion and group delay properly. next is to implement more complexes acoustica loads.
here is a correct code for those who are curious :
Code:
% General constants
c=343.29; rho=1.2039; P_ref=20*(10^(-6)); Air_Imp=413.27;
% frequencies and wave number
f=logspace(1,4,300); w=2*pi*f; %frequencies [Hz]
wc = w / c; % wave number
% distance and input voltage
r = 1; % Distance to from the source [m]
eg = 2.83; % Input voltage [V]
% Thiele/Small parameters for RCF MB15N405
Sd=0.09; fs=46; Vas=0.124; Qms=4.8; Bl=23.5; Re=5.5; Le=0.0011; Mms=0.098; Rms=6.74;
Xmax=7; Qts=0.27; Qes=0.28;
% Radiation impedance uing bessel and struve functions choose
[zar,zmr]=RadImpC(Sd); % zar for acoustic domain and zmr for mechanical domain
za=2*zar; % takes into account both sides of the driver
% Calculating mechanical mass Mmd
Mr=imag(zmr)/(w); % Radiation mass calculated from Kolbrek eq.14.4.5
Mmd=Mms-Mr; % Moving mass without radiation mass [kg]
% Calculating mechanical compliance Cms
Cms = (Vas/(rho*(c^2)))/(Sd^2); % Cms from Kolbrek eq.16.7.6 [m/N]
% acoustical equivalence
Pg=eg*Bl/(Sd*Re); % equivalent pressure [Pa]
mas=Mms/(Sd^2); % equivalent mass [kg/m^4]
% Loop for each frequency
for k = 1:length(w)
% Transfer matrices for sub domains
tre(:, :, k) = [1, Re; 0, 1]; % voice coil resistance
tle(:, :, k) = [1, 1i*w(k)*Le; 0, 1]; % Inductance
tem(:, :, k) = [0, Bl; 1/Bl, 0]; % gyrator
tmm(:, :, k) = [1, 1i*w(k)*Mmd; 0, 1]; % Mobil mass
tcm(:, :, k) = [1, 1/(1i*w(k)*Cms); 0, 1]; % Compliance
trm(:, :, k) = [1, Rms; 0, 1]; % mechanic losses
% Acoustic domain
tma(:, :, k) = [Sd, 0; 0, 1/Sd]; % Tranformer
tza(:, :, k) = [1, za(k); 0, 1]; % Radiation Impedance
% Complete chain from eletrical domain to acoustical domain
Tt(:, :, k) = tre(:, :, k) * tle(:, :, k) * tem(:, :, k) * tmm(:, :, k) * tcm(:, :, k) * trm(:, :, k) * tma(:, :, k) * tza(:, :, k);
% Calculating Iin [A] electrical current from the source from
% electro-mechanico-acoustical matrix
Pu_forced = [1; 1]; % imposes a specific output to force calcul of input current
In = Tt(:,:,k) * Pu_forced; % calculating inputs vector
zin = In(1) / In(2); % calculating impedance seen from the source
iin=eg/zin; % calculating input current from zin
In = [eg;iin]; % electrical nput vector
% Calculating the final response
Ttinv = inv(Tt(:, :, k)); % Inverting matrix
Pu = Ttinv * In; % Calculating results
pout(k) = Pu(1); % Acoustic pressure at Sd
uout(k) = Pu(2); % Acoustic volume velocity
endSPL, phase, excursion and group delay can be calculated from uout easily.
some curves bellow, in red Hornresp and in blue my model, we can observe a perfect matching !
thanks to everyone who took time reading this thread.
Last edited: