%   Script LAUTSPR.M zur Herstellung von Bildern
%
%   Stabilitätsanalyse der Rückkopplung vom Lautsprecher
%   zum Mikrofon
%
% J. Lunze
% 17.11.1995
% für 2. Auflage: 7.12.1998
% für 3. Auflage: 3.4.2001
% für 4. Auflage: 18.9.2003
% für 5. Auflage: 18.4.2005
% für 6. Auflage: 4.12.2006
% für 7. Auflage: 28.4.2008
% für 9. Auflage: 13.2.2012
% für 10. Auflage: 3.2.2014
% für 12. Auflage: 5.08.2019
%
%   
clear
close all

k_V=1000;        %             Verstärker
wgr_V = 125000;  % rad/s
k_L = 1;         % N/m^2/V     Lautsprecher
wgr_L = 60000;   % rad/s
%d=2;            % m           Entfernung
c=335;           % m/s         Schallgeschwindigkeit
%Tt = d/c;       % s              -> Totzeit
%k_W=1/d;        %                   Abnahme des Schalldrucks
k_M = 0.02;      % V m^2/N     Empfindlichkeit des Mikrofons
wgr_M = 60000;   % rad/s
k_W=5*5/2;       %             Proportionalitätsfaktor für
                 %             Abnahme des Schalldrucks

%%    Frequenzgang des offenen Lautsprecherkreises
%
figure(1);
subplot(3,3,[1,2,4,5]);
Pole=[-wgr_V -wgr_L -wgr_M];
n=poly(Pole)/wgr_V/wgr_L/wgr_M;
for d=[10, 15, 25]
   Tt = d/c;
   k_W1=k_W/d/d;
   z=[k_L*k_W1*k_M*k_V];
   System=tf(z, n);
   % Frequenz für Phasenverschiebung 360 Grad
   w360=2*pi/Tt;
   W = logspace(0, log10(w360), 600)';
   [Betrag, Phase] = bode(System, W);
   for k=1:length(W)
      Phase(k) = Phase(k)-180 - W(k)*Tt/pi*180;
        % Minuszeichen vor G0 durch Phasenverschiebung
        % von -180 Grad berücksichtigt
      Real(k) = Betrag(k)*cos(Phase(k)*pi/180);
      Imag(k) = Betrag(k)*sin(Phase(k)*pi/180);
   end
   plot(Real, Imag, 'b');
   hold on;
end;
plot([-1 -1],[-0.2 0.2],'k',[-1.2 -0.8],[0 0],'k');
latextitle('Ortskurve des Mikrofon-Lautsprecher-Systems');
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis('square');
axis([-3, 3, -3, 3]);
xticks([-2 0 2]);
yticks([-2 0 2]);
grid('on');
latextext(0.1, 2.7, '$$d$$=10m');
latextext(0.5, 0.4, '25m');
latextext(0.5, 1.4, '15m ');
hold off



%%     Amplitudengang der offenen Kette für d=15
%
figure(2)
W = logspace(0, 7, 800)';
d=15;
Tt = d/c;
k_W1=k_W/d/d;
z=[k_L*k_W1*k_M*k_V];
System=tf(z, n);
[Betrag15, Phase15] = bode(System, W);
for k=1:size(W,1)
   Phase15(k) = Phase15(k) -180 - W(k)*Tt/2/pi*360;
end;
subplot(2,1,1)
semilogx(W, 20*log10(Betrag15(:)),'b');
latextitle('Lautsprecher-Mikrofon-Anordnung $$(d=15m)$$');
xleer;
latexylabel('$$|G|$$ in dB');
rotateY;
grid('on');
axis([1 10000000 -100 20]);
yticks([-100 -50 0]);
%
subplot(2,1,2);
semilogx(W, Phase15(:),'b');
latexylabel('$$\phi$$ in Grad');
rotateY;
grid('on');
latexxlabel('$$\omega$$ in rad/s');
axis([1 10000000 -450 0]);
xticks([1 100 10000 1000000 100000000]);
yticks([-360 -180 0]);
 

%%    Analyse des Rückkopplungskreises beim Telefon
%
figure(3)
subplot(3,3,[1,2,4,5]);
k_Vk_L = 10;      % N/m^2/V     Lautsprecher
d=0.3;            % m           Entfernung
Tt = d/c;         % s              -> Totzeit
k_W=1*1/2;        %             Proportionalitätsfaktor
k_W1=k_W/d/d;     %                   Abnahme des Schalldrucks
wgr_M = 4000;     % rad/s
%   Frequenzgang des offenen Lautsprecherkreises   
Pole=[-wgr_V -wgr_L -wgr_M];
n=poly(Pole)/wgr_V/wgr_L/wgr_M;
z=[k_Vk_L*k_W1*k_M];
System=tf(z, n);
% Frequenz für Phasenverschiebung 360 Grad
w360=2*pi/Tt;
W = logspace(0, log10(2*w360), 600)';
[Betrag, Phase] = bode(System, W);
Phase(:) = Phase(:) -180 - Tt/pi*180*W(:);
Real = Betrag(:).*cos(Phase(:)*pi/180);
Imag = Betrag(:).*sin(Phase(:)*pi/180);
plot(Real(:), Imag(:), 'b');
hold on;
plot([-1 -1],[-0.1 0.1],'k', [-1.1 -0.9],[0 0],'k');
latextitle('Ortskurve des Telefon-Lautsprecher-Systems');
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis('square');
axis([-1.2, 1.2, -1.2, 1.2]);
xticks([-1 0 1]);
yticks([-1 0 1]);
grid('on');
latextext(1.02, 0.1, '$$\downarrow$$');
hold off


%%     Bodediagramm für diese Anordnung
%
figure(4)
subplot(2,1,1);
semilogx(W, 20*log10(Betrag(:)),'b');
xleer;
latexylabel('$$|G|$$ in dB');
rotateY;
grid('on');
latextitle('Telefon-Lautsprecher-System');
axis([10 10000 -20 10]);
yticks([-20 -10 0 10]);
%
subplot(2,1,2);
semilogx(W, Phase(:),'b');
latexxlabel('$$\omega$$ in rad/s');
latexylabel('$$\phi$$ in Grad');
rotateY;
grid('on');
axis([10 10000 -360 0]);
yticks([-270 -180 -90 0]);
xticks([10 100 1000 10000]);

%%     Bilder auf .EPS- und .PS-Dateien
%
Epsbildc15('lautspr');
