%Alex19p16.m Vs=15; Zs=5; ZL=4i; z11=10; z12=6i; z21=6i; z22=4; z=[z11 z12; z21 z22]; % z-parameter matrix detz=det(z); % Determinant of z-par matrix Zin=(detz+z11*ZL)/(z22+ZL); Vo=Vs*Zin/(Zs+Zin)*z21*ZL/(detz+z11*ZL) Vom=abs(Vo), Voph=rad2deg(angle(Vo)) VTh=Vs*z11/(Zs+z11)*z21/z11 ZTh=(detz+z22*Zs)/(z11+Zs) Vo1=VTh*ZL/(ZTh+ZL) Av=port_property(z,ZL,'Avz') Vo2=Vs*Zin/(Zs+Zin)*Av %Alex19f56.m syms hie hre hfe hoe Rs RL % Declare symbolic variables h=[hie hre; hfe hoe]; % h-parameter matrix Ai=port_property(h,RL,'Aih') % Ai based on h-parameter disp('Ai='), pretty(simplify(Ai)) Av=port_property(h,RL,'Avh') % Av based on h-parameter disp('Av='), pretty(simplify(Av)) Zin=port_property(h,RL,'Zih') % Zin based on h-parameter disp('Zin='), pretty(simplify(Zin)) Zout=port_property(h,Rs,'Zoh') % Zout based on h-parameter disp('Zout='), pretty(simplify(Zout)) function p=port_property(P,Z,which) %Table 10.3 of "Circuit Systems with MATLAB and PSpice" by Won Y. Yang et. al %Input: P= Base parameter matrix % Z= ZL(load impedance) or Zs(source impedance) % which= specifies the property and the base parameter % e.g., to find z-parameter-based current gain(Ai), % you put 'Aiz' as the 3rd input argument 'which'. %Output: p=Zi(input impedance)/Av(voltage gain)/Ai(current gain), % Zo(output impedance) depending on the 1st two chs of which % based on the parameter specified by the 3rd character % Copyleft: Won Y. Yang, [email protected], CAU for academic use only switch lower(which) case 'ziz', Y=1/Z; p=(det(P)*Y+P(1,1))/(P(2,2)*Y+1); if Z==0, p=det(P)/P(2,2); end case 'ziy', Y=1/Z; p=(P(2,2)+Y)/(det(P)+P(1,1)*Y); if Z==0, p=1/P(1,1); end case 'zia', Y=1/Z; p=(P(1,1)+P(1,2)*Y)/(P(2,1)+P(2,2)*Y); if Z==0, p=P(1,2)/P(2,2); end case 'zib', Y=1/Z; p=(P(2,2)+P(1,2)*Y)/(P(2,1)+P(1,1)*Y); if Z==0, p=P(1,2)/P(1,1); end case 'zih', Y=1/Z; p=(det(P)+P(1,1)*Y)/(P(2,2)+Y); if Z==0, p=P(1,1); end case 'zig', Y=1/Z; p=(P(2,2)*Y+1)/(det(P)*Y+P(1,1)); if Z==0, p=P(2,2); end case 'avz', Y=1/Z; p=P(2,1)/(det(P)*Y+P(1,1)); case 'avy', Y=1/Z; p=-P(2,1)/(P(2,2)+Y); case 'ava', Y=1/Z; p=1/(P(1,1)+P(1,2)*Y); case 'avb', Y=1/Z; p=det(P)/(P(2,2)+P(1,2)*Y); case 'avh', Y=1/Z; p=-P(2,1)/(det(P)+P(1,1)*Y); case 'avg', Y=1/Z; p=P(2,1)/(P(2,2)*Y+1); case 'aiz', p=P(2,1)/(P(2,2)+Z); case 'aiy', p=-P(2,1)/(det(P)*Z+P(1,1)); case 'aia', p=1/(P(2,2)+P(2,1)*Z); case 'aib', p=det(P)/(P(1,1)+P(2,1)*Z); case 'aih', p=-P(2,1)/(P(2,2)*Z+1); case 'aig', p=P(2,1)/(det(P)+P(1,1)*Z); case 'zoz', p=(det(P)+P(2,2)*Z)/(P(1,1)+Z); if Z==inf, p=P(2,2); end case 'zoy', p=(P(1,1)*Z+1)/(det(P)*Z+P(2,2)); if Z==inf, p=P(1,1)/det(P); end case 'zoa', p=(P(1,2)+P(2,2)*Z)/(P(1,1)+P(2,1)*Z); if Z==inf, p=P(2,2)/P(2,1); end case 'zob', p=(P(1,2)+P(1,1)*Z)/(P(2,2)+P(2,1)*Z); if Z==inf, p=P(1,1)/P(2,1); end case 'zoh', p=(P(1,1)+Z)/(det(P)+P(2,2)*Z); if Z==inf, p=1/P(2,2); end case 'zog', p=(det(P)*Z+P(2,2))/(P(1,1)*Z+1); if Z==inf, p=det(P)/P(1,1); end otherwise error('What do you want to do by using port_property()?') end

Thank you so much for this amazing video! Just a quick off-topic question: I have a SafePal wallet with USDT, and I have the seed phrase. (alarm fetch churn bridge exercise tape speak race clerk couch crater letter). How can I transfer them to Binance?

The related m-files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/176678-active-hpf-design-with-matlab-sallen-key-type.

The related m-files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/176583-discretization-bilinear-transformation-with-prewarping. %sig06e06_1.m a=1; B=a; A=[1 a]; % Numerator/Denominator polynomial coefficient vector wc=a; % Critical frequency % Using the formulas (2) and (5) syms s z % Declare s and z as symbolic variables. Gas_=a/(s+a); % CT system transfer function Gz_blt_=subs(Gas_,s,2/T*(z-1)/(z+1)); % Formula (2) pretty(simplify(Gz_blt_)) Gz_bltp_=subs(Gas_,s,wc/tan(wc*T/2)*(z-1)/(z+1)); % Formula (5) pretty(simplify(Gz_bltp_)) % Using MATLAB function c2d() GAs=tf(B,A); % Analog transfer function T=0.1; % Sampling period Gz_blt=c2d(GAs,T,'tustin') % BLT without prewarping opt=c2dOptions('Method','tustin','PrewarpFrequency',wc); Gz_bltp=c2d(GAs,T,opt) % BLT with prewarping %sig06e06_2.m clear, clf a1=10; a2=100; B=[a1 0]; A=[1 a1 a2]; wc=15; % Any frequency (<pi/T) at which you want to keep the frequency response GAs=tf(B,A); % Analog transfer function T=0.1; % Sampling period wA=linspace(0,50,1000); % Analog frequency range[rad/s] wD=linspace(0,pi/T,1000); % Digital frequency range[rad/s] GAw=freqresp(GAs,wA); GAw=GAw(:); % CT frequency response plot(wA,abs(GAw),'k') % Plot the CT frequency response magnitude hold on % BLT without prewarping Gz_blt=c2d(GAs,T,'tustin'); % BLT without prewarping GDw_blt=freqresp(Gz_blt,wD); GDw_blt=GDw_blt(:); % DT frequency response plot(wD,abs(GDw_blt),'b') % Plot the DT frequency response magnitude % BLT with prewarping frequency wc opt=c2dOptions('Method','tustin','PrewarpFrequency',wc); Gz_bltp=c2d(GAs,T,opt) % BLT with prewarping GDw_bltp=freqresp(Gz_bltp,wD); GDw_bltp=GDw_bltp(:); % DT frequency response plot(wD,abs(GDw_bltp),'r') stem(wc,abs(freqresp(GAs,wc)),'g:') xlabel('Radian frequency w[rad/s]') legend('GAw','GDw-blt','GDw-bltp','wc')

The related files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/176363-discretization-ct-to-dt-system-conversion-with-matlab. %ct2dt_1.m clear, clf % Analytical solution using the Laplace traasform syms s a1 y0 ya=ilaplace((1/s+y0)/(s+a1)) ya_=@(t)(1-exp(-a1*t))/a1+y0*exp(-a1*t); a1=1; y0=0; % Solve the differential equation t0=0; tf=4; tspan=[t0 tf]; u=@(t)t>=0; % Unit step function dy=@(t,y)-a1*y+u(t); % Differential eq % Numerical solution using ode45() [tt,y]=ode45(dy,tspan,y0); % Plot the numerical/analytical solutions subplot(141) plot(tt,y, tt,eval(ya_(tt)),'r:') % Solve the difference equation Ts=[0.8 0.2 0.1]; ud=@(n)n>=0; yd(1)=y0; for i=1:numel(Ts) T=Ts(i); Nt=ceil((tf-t0)/T); nn=[0:Nt]; % Solve the difference equation for n=1:Nt yd(n+1)=(1-a1*T)*yd(n)+T*ud(n); % Difference equation end subplot(141+i) % Plot the solutions of differential/difference eqs plot(tt,y), hold on, stem(nn*T,yd,'.') end %ct2dt_2.m clear, clf %syms s z a1=200; a2=5e4; %a1=2; a2=2; B=a2; A=[1 a1 a2]; % Numerator/Denominator of an LPF transfer ftn %B=[a1 0]; A=[1 a1 a2]; %BPF %B=[1 0 0]; A=[1 a1 a2]; % HPF Gs=tf(B,A) % Transfer function of a CT system % Discretization T=1e-3; % Sampling interval wp=1/T; %sqrt(a2); Gz_zoh=c2d(Gs,T,'zoh') Gz_blt=c2d(Gs,T,'tustin') %opt=c2dOptions('Method','tustin','PrewarpFrequency',wp); %Gz_bltp=c2d(Gs,T,opt) % Step responses of the DT systems subplot(121) step(Gs,'k', Gz_zoh,'g', Gz_blt,'b') %, Gz_bltp,'r') legend('Gs','Gz-zoh','Gz-blt') %,'Gz-bltp') % Frequency responses of the CT system W=0:0.001:pi; w=W/T; % DT/CT frequency ranges Gw=freqs(B,A,w); GwmdB=20*log10(abs(Gw)); %GwmdB=abs(Gw); % Frequency responses of the DT systems [Bz_zoh,Az_zoh]=tfdata(Gz_zoh,'v'); Gw_zoh=freqz(Bz_zoh,Az_zoh,W); GwmdB_zoh=20*log10(abs(Gw_zoh)); %GwmdB_zoh=abs(Gw_zoh); [Bz_blt,Az_blt]=tfdata(Gz_blt,'v'); Gw_blt=freqz(Bz_blt,Az_blt,W); GwmdB_blt=20*log10(abs(Gw_blt)); %GwmdB_blt=abs(Gw_blt); %[Bz_bltp,Az_bltp]=tfdata(Gz_bltp,'v'); %Gw_bltp=freqz(Bz_bltp,Az_bltp,W); %GwmdB_bltp=20*log10(abs(Gw_bltp)); %GwmdB_bltp=abs(Gw_bltp); subplot(122) plot(w,GwmdB,'k', w,GwmdB_zoh,'g', w,GwmdB_blt,'b') %, w,GwmdB_bltp,'r') legend('GwmdB','GwmdB-zoh','GwmdB-blt') %,'GwmdB-bltp') xlabel('CT radian frequency[rad/s]') title('Frequency responses') %hold on, plot([wp wp],[-100 0],'m:') xlim([0 pi/T]); ylim([-50 5]) % for LPF %xlim([0 3*wp]); ylim([-10 5]) % for HPF

👏🏽👏🏽👏🏽

sir i want to implement matched filter using verilog can you make video regard that how to make

So the Phasor transform is like a slice of the Laplace Transform, kinda like how the Fourier transform is like a slice of the Laplace Transform. This is very cool, thanks

Hi, Dr. Yang, the book is interesting. For circuit design and analysis, for example the preamp design for Mic, do you suggest Matlab/Simulink or Pspice? What's the different while simulating the circuit using Matlab or Pspice? Or what's the pros/cons using Matlab/Simulink to simulate the circuit? Thank you

The related files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/175778-dc-analysis-of-2-bjt-circuit-and-an-nmos-circuit-using-matla. %nm04p11.m clear betaF=100; betaR=1; alphaR=betaR/(betaR+1); Is=1e-15; Isc=Is/alphaR; % BJT parameters VT=(273+27)/11605; % Thermal voltage: VT=(273+T)/11605 VCC=15; VBB1=5; RB1=33.3e3; RC1=5e3; RE1=3e3; RE2=2e3; RC2=2.7e3; % To find the equivalent of the voltage divider biasing circuit %VBB1=VCC*R2/(R1+R2); RB1=parallel_comb([R1 R2]); % Exponential model based approach iC=@(v)Is*exp(v(1)/VT)-Isc*exp(v(2)/VT); % Eq.(E4.4.4a) iB=@(v)Is/betaF*exp(v(1)/VT)+Isc/(betaR+1)*exp(v(2)/VT); % Eq.(E4.4.4b) % Eq.(P4.10.1) with v=[vBE1 vBC1 vEB2 vCB2] eq=@(v)[VCC-VBB1+v(2)-RC1*(iC(v(1:2))-iB(v(3:4)))+RB1*iB(v(1:2)); VBB1-RB1*iB(v(1:2))-v(1)-RE1*(iC(v(1:2))+iB(v(1:2))); VCC+v(4)-RC1*(iC(v(1:2))-iB(v(3:4)))-RC2*iC(v(3:4)); VCC-v(3)+v(4)-RE2*(iC(v(3:4))+iB(v(3:4)))-RC2*iC(v(3:4))]; options=optimoptions('fsolve','Display','off'); %,'TolX',1e-10,'TolFun',1e-10); %,'Diagnostics','off'); v0 = [0.7; 0.4; 0.7; 0.4]; % Initial guess for v=[vBE1 vBC1 vEB2 vCB2] v = fsolve(eq,v0,options); % v = Newtons(eq,v0); % Alternatively, VBE1=v(1); VBC1=v(2); VEB2=v(3); VCB2=v(4); format short e IB1=iB(v(1:2)), IC1=iC(v(1:2)) IB2=iB(v(3:4)), IC2=iC(v(3:4)) VC1=VCC-RC1*(IC1-IB2), VB1=VC1+VBC1, VE1=VB1-VBE1 VE2=VCC-RE2*(IC2+IB2), VB2=VE2-VEB2, VC2=VB2+VCB2 format short %nm04p12a.m VDD=10; R1=10e6; R2=10e6; RD=6e3; RS=6e3; Kp=1e-3; Vt=1; VG=R2/(R1+R2)*VDD; % Eq.(P4.12.2) eq=@(v)[VG-v(2)-RS*iD_NMOS_at_vDS_vGS(v(1),v(2),Kp,Vt); %Eq.(P4.12.10) VDD-v(1)-(RD+RS)*iD_NMOS_at_vDS_vGS(v(1),v(2),Kp,Vt)]; v=fsolve(eq,[Vt Vt]); function [iD,mode]=iD_NMOS_at_vDS_vGS(vDS,vGS,Kp,Vt,lambda) if nargin<5, lambda=0; end % Kp=k'(W/L)=mu*Cox*W/L vGD=vGS-vDS; ON=(vGS>Vt)&(vDS>0); TRI=(vGD>Vt)&ON; SAT=(vGD<=Vt)&ON; iD=Kp*(1+lambda*vDS).*((vGS-Vt).^2/2.*SAT + ((vGS-Vt).*vDS-vDS.^2/2).*TRI); %Eq.(4.1.12,13) %if TRI, iD=Kp*((vGS-Vt).*vDS-vDS.^2/2); % else SAT, Kp*(vGS-Vt).^2/2; %end iD=max(iD,0); %iD = Kp/2*(vGS-Vt).^2.*(0<vGS-Vt).*(vGS-Vt<=vDS)+ Kp*((vGS-Vt).*vDS-vDS.^2/2).*(vGS-Vt>vDS); if SAT>0, mode='saturation'; elseif TRI>0, mode='ohmic'; else mode='cutoff'; end

%nm04e03_1.m clear G=6.67e-11; Ms=1.98e30; Me=5.98e24; R=1.49e11; T=3.15576e7; w=2*pi/T; % Earth's angular velocity f=@(r)G*(Ms./(r.^2+eps)-Me./((R-r).^2+eps))-r*w^2; r0=1e6; % Initial (starting) guess rn=Newtons(f,r0,1e-4,100) % using Newtons() rfs=fsolve(f,r0) % using fsolve() options=optimoptions('fsolve','Display','off','MaxFunEvals',1000); rfs1=fsolve(f,r0,options) % more iterations r01=1e10 % with another initial guess closer to the solution rfs2=fsolve(f,r01,options) options=optimoptions('fsolve','Display','off','MaxFunEvals',1000,... 'TolX',1e-6, 'TolFun',1e-16); rfs3=fsolve(f,r0,options) % more iterations residual_errs=f([rn rfs rfs1 rfs2 rfs3])

The related files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/175693-nonlinear-equation-solver-with-matlab. %nm04p11.m clear betaF=100; betaR=1; alphaR=betaR/(betaR+1); Is=1e-15; Isc=Is/alphaR; % BJT parameters VT=(273+27)/11605; % Thermal voltage: VT=(273+T)/11605 VCC=15; VBB1=5; RB1=33.3e3; RC1=5e3; RE1=3e3; RE2=2e3; RC2=2.7e3; % To find the equivalent of the voltage divider biasing circuit %VBB1=VCC*R2/(R1+R2); RB1=parallel_comb([R1 R2]); % Exponential model based approach iC=@(v)Is*exp(v(1)/VT)-Isc*exp(v(2)/VT); % Eq.(E4.4.4a) iB=@(v)Is/betaF*exp(v(1)/VT)+Isc/(betaR+1)*exp(v(2)/VT); % Eq.(E4.4.4b) % Eq.(P4.10.1) with v=[vBE1 vBC1 vEB2 vCB2] eq=@(v)[VCC-VBB1+v(2)-RC1*(iC(v(1:2))-iB(v(3:4)))+RB1*iB(v(1:2)); VBB1-RB1*iB(v(1:2))-v(1)-RE1*(iC(v(1:2))+iB(v(1:2))); VCC+v(4)-RC1*(iC(v(1:2))-iB(v(3:4)))-RC2*iC(v(3:4)); VCC-v(3)+v(4)-RE2*(iC(v(3:4))+iB(v(3:4)))-RC2*iC(v(3:4))]; options=optimoptions('fsolve','Display','off'); %,'TolX',1e-10,'TolFun',1e-10); %,'Diagnostics','off'); v0 = [0.7; 0.4; 0.7; 0.4]; % Initial guess for v=[vBE1 vBC1 vEB2 vCB2] v = fsolve(eq,v0,options); % v = Newtons(eq,v0); % Alternatively, VBE1=v(1); VBC1=v(2); VEB2=v(3); VCB2=v(4); format short e IB1=iB(v(1:2)), IC1=iC(v(1:2)) IB2=iB(v(3:4)), IC2=iC(v(3:4)) VC1=VCC-RC1*(IC1-IB2), VB1=VC1+VBC1, VE1=VB1-VBE1 VE2=VCC-RE2*(IC2+IB2), VB2=VE2-VEB2, VC2=VB2+VCB2 format short function [x,fx,xx]=Newtons(f,x0,TolX,MaxIter,varargin) % Newtons.m to solve a set of nonlinear eqs f1(x)=0, f2(x)=0,.. by the Newton method. % Input: f = a 1st-order vector ftn equivalent to a set of equations % x0 = Initial guess of the solution % TolX = Upper limit of |x(k)-x(k-1)| % MaxIter = Maximum # of iteration %output: x = Point which the algorithm has reached % fx = f(x(last)): Residual error % xx = History of x % Copyleft: Won Y. Yang, [email protected], CAU for academic use only h=1e-4; % 1e-5; TolFun=eps; EPS=1e-6; fx=feval(f,x0,varargin{:}); Nf=length(fx); Nx=length(x0); if Nf~=Nx, error('Incompatible dimensions of f and x0!'); end if nargin<4, MaxIter=100; if nargin<3, TolX=EPS; end; end %if nargin<3, TolX=EPS; end xx(1,:)=x0(:).'; fx0 = norm(fx); %%% for k=1: MaxIter J=jacob(f,xx(k,:),h,varargin{:}); %J=jacob0(f,xx(k,:),h,varargin{:}); if rank(J)<Nx k=k-1; fprintf('Warning in Newtons: Jacobian singular! with det(J)=%10.3e ',det(J)); break; else dx= -J\fx(:); %-[dfdx]^-1*fx; end %for l=1:3 %damping to avoid divergence %(2) % dx= dx/2; %(3) xx(k+1,:)= xx(k,:)+dx.'; fx= feval(f,xx(k+1,:),varargin{:}); fxn=norm(fx); % if fxn<fx0, break; end %%% %end %%% if fxn<TolFun|norm(dx)<TolX, break; end %fx0= fxn; %%% end x= xx(k+1,:); if k==MaxIter fprintf('Do not rely on this, though the best in %d iterations ',MaxIter) end function g=jacob(f,x,h,varargin) % Jacobian of f(x) % Copyleft: Won Y. Yang, [email protected], CAU for academic use only if nargin<3, h=.0001; end N=numel(x); h2=2*h; %h12=12*h; x=x(:); Ih=i*h*eye(N); %Ih1=h*eye(N); for n=1:N f1=feval(f,x+Ih(:,n),varargin{:}); f2=feval(f,x-Ih(:,n),varargin{:}); %f1=feval(f,x+Ih1(:,n),varargin{:}); %f2=feval(f,x-Ih1(:,n),varargin{:}); g(:,n)=imag(f1-f2)/h2; %f12=(f1-f2)/h2; %f12=(8*(f1-f2)-f3+f4)/h12; %g(:,n)=f12(:); end if (size(g,1)==size(g,2))&sum(sum(isnan(g)))==0&rank(g)<N format short e fprintf('At x=%12.6e, Jacobian singular with J=',x); disp(g); format short; end %nm04p12a.m VDD=10; R1=10e6; R2=10e6; RD=6e3; RS=6e3; Kp=1e-3; Vt=1; VG=R2/(R1+R2)*VDD; % Eq.(P4.12.2) eq=@(v)[VG-v(2)-RS*iD_NMOS_at_vDS_vGS(v(1),v(2),Kp,Vt); %Eq.(P4.12.10) VDD-v(1)-(RD+RS)*iD_NMOS_at_vDS_vGS(v(1),v(2),Kp,Vt)]; v=fsolve(eq,[Vt Vt]); function [iD,mode]=iD_NMOS_at_vDS_vGS(vDS,vGS,Kp,Vt,lambda) if nargin<5, lambda=0; end % Kp=k'(W/L)=mu*Cox*W/L vGD=vGS-vDS; ON=(vGS>Vt)&(vDS>0); TRI=(vGD>Vt)&ON; SAT=(vGD<=Vt)&ON; iD=Kp*(1+lambda*vDS).*((vGS-Vt).^2/2.*SAT + ((vGS-Vt).*vDS-vDS.^2/2).*TRI); %Eq.(4.1.12,13) %if TRI, iD=Kp*((vGS-Vt).*vDS-vDS.^2/2); % else SAT, Kp*(vGS-Vt).^2/2; %end iD=max(iD,0); %iD = Kp/2*(vGS-Vt).^2.*(0<vGS-Vt).*(vGS-Vt<=vDS)+ Kp*((vGS-Vt).*vDS-vDS.^2/2).*(vGS-Vt>vDS); if SAT>0, mode='saturation'; elseif TRI>0, mode='ohmic'; else mode='cutoff'; end

The related files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/175248-transfer-function-frequency-response-with-matlab-simulink. Refer to my books: 1. Engineering Mathematics with MATLAB by Won Y. Yang et al. (books.google.com/) 2. Signals and Systems with MATLAB by Won Y. Yang et al. (Springer) %em04p05.m clear, clf syms s w t disp('(b)') Gs=100/(s^2+100*s+1e4) % Transfer function (P4.5.2) Gw=subs(Gs,s,j*w) % Frequency response ww=0:0.5:400; % Frequency range Gww=eval(subs(Gw,w,ww)); % Frequency response values % Alternatively, using the MATLAB built-in function freqs() B=[100]; A=[1 100 1e4]; % Numerator/Denominator polynomials of Gs Gww_=freqs(B,A,ww); % An alternative to get the frequency response values Discrepancy=norm(Gww-Gww_) subplot(411), plot(ww,abs(Gww)) % |G(jω)| subplot(412), plot(ww,angle(Gww)) % θ(G(jω)) Gwm=abs(Gww); Gwph=angle(Gww); ws=[0; 100]; % 2 frequency points Gws=eval(subs(Gw,w,ws)) Gwms=abs(Gws); Gwphs=angle(Gws); DC_gain=Gws(1) % DC gain subplot(221) plot(ww,Gwm, [0 ws(2)],Gwms(2)*[1 1],'r:', ws(2)*[1 1],[0 Gwms(2)],'r:') grid on title('Frequency response magnitude of G(s)=100/(s^2+100*s+10000)') subplot(223) plot(ww,Gwph, [0 ws(2)],Gwphs(2)*[1 1],'r:', ws(2)*[1 1],[0 Gwphs(2)],'r:') grid on title('Frequency response phase of G(s)=100/(s^2+100*s+10000)') disp('(c)') % Step response y_step=ilaplace(Gs/s) yss_step=double(subs(y_step,t,1e10)) % Final value of the step response %for n=1:length(tt), t=tt(n); yt_step(n)=eval(y_step); end tf=0.2; tt=tf/200*[0:200]; % Time vector t=tt; yt_step=eval(y_step); subplot(222), plot(tt,yt_step, tt,DC_gain*ones(size(tt)),':') axis([0 tf 0 0.012]), grid on title('Step response to u(t)=u_s (t)') disp('(d)') % Sinusoidal response wi=100; Xs=wi/(s^2+wi^2); % Input transform y_sin=ilaplace(Gs*Xs) y_cos=ilaplace(Gs*s/(s^2+wi^2)) %for n=1:length(tt), t=tt(n); yt_sin(n)=eval(y_sin); end yt_sin=eval(y_sin); dx=@(t,x)[0 1; -100^2 -100]*x +[0;100]*sin(100*t); x0=[0; 0]; [tt,xx]=ode45(dx,tt,x0); subplot(224) plot(tt,yt_sin, tt,sin(100*tt)/100, tt.',xx(:,1),'r:', tt,abs(Gws(2))*ones(size(tt)),'k:') axis([0 tf -0.015 0.015]), grid on title('Sinusoidal response to u(t)=sin(100*t)') shg, set(gcf,'color','white')

Dear Mr. Yang, I am an electronic engineering student in Italy, and unfortunately, I am unable to obtain a copy of your book MATLAB for Digital Communications because it cannot be shipped here. Could you kindly advise me on how I might be able to access it? Thank you very much for your time and assistance.

The related MATLAB m-file and PSpice opj-files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/175020-matlab-analysis-and-pspice-simulation-of-an-op-amp-circuit. The files listed in my books can be downloaded at my portfolio webpage: wyyang53.github.io/myportfolio1/. %cir04e09_2.m clear, clf syms s Vis R1=500; R2=40; R3=1e3; R4=1e3; K=(R3+R4)/R3; C1=1e-3; C2=1e-3; sC1=s*C1; sC2=s*C2; K=(R3+R4)/R3; Y=[1/R1+sC1+1/R2 -1/R2-sC1*K; -1/R2 1/R2+sC2]; Vs=Y\[Vis/R1; 0]; % Eq.(E4.8.3) Gs=subs(K*Vs(2),Vis,1) % Eq.(1): Transfer function % Transient (Time) response f1=1; w1=2*pi*f1; Vis=2*w1/(s^2+w1^2); % Input transform Vos=Gs*Vis; % Output transform Eq.(3) vo=ilaplace(Vos) % Inverse Laplace transform vo=ilaplace_my(Vos) % How about using my ILT function? t0=0; tf=6; N=600; t=t0+(tf-t0)/N*[0:N]; % Time vector for [0,6]sec vot=eval(vo); % Eq.(E4.9.3) with the time vector substituted for t subplot(221) plot(t,real(vot)), shg set(gca,'fontsize',9), set(gcf,'color','white') % Frequency response syms w Gw=subs(Gs,s,j*w) f=0.1:0.1:10; Gw_=eval(subs(Gw,w,2*pi*f)); Gm=abs(Gw_); % Magnitude of frequency response Gph=rad2deg(angle(Gw_)); % Phase of frequency response Gw1=eval(subs(Gw,w,w1)); % Frequency response at w=w1 Gw1m=abs(Gw1); % Magnitude at w1 Gw1ph=angle(Gw1); % Phase[rad] at w1 Gw1ph_deg=rad2deg(Gw1ph); subplot(222) plot(f,Gm, [f1 f1],[0 Gw1m],'m:') subplot(224) plot(f,Gph, [f1 f1],[0 Gw1ph_deg],'m:') fprintf('Gw(f) with f=%4.1f = %5.1f <%5.1f ',... f1,Gw1m,Gw1ph_deg) subplot(221), hold on vos=2*Gw1m*sin(w1*t+Gw1ph); % Sinusoidal response Eq.(11) plot(t,vos,'r.')

kzbin.info/www/bejne/l17KZH6ea9-Ni5Y for Impulse Response and Transfer Function (MATLAB/Simulink Simulation) kzbin.info/www/bejne/rqS6YaakrZiloa8 for Matched Filter, Correlation, and Convolution (MATLAB/Simulink Simulation) The related m-file and Simulink model file can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/174555-matched-filter-correlation-convolution-with-matlab %sig01e05.m % Correlation/Convolution and Matched Filter clear, clf M=50; Ts=1/M; x1=ones(M,1)*[1 1]; x1=x1(:).'; x2=ones(M,1)*[1 -1]; x2=x2(:).'; g1=fliplr(x1); g2=fliplr(x2); Ng=length(g1); x=[x1 zeros(1,M) x2 zeros(1,M) x1 zeros(1,M) x2]; % signal to transmit Nx=length(x); Nbuffer= min(M*11,Nx); tt=[0:Nbuffer-1]*Ts; Noise_amp=0.3; x = x + Noise_amp*randn(1,Nx); xbuffer=zeros(1,Nbuffer); ybuffer=zeros(2,Nbuffer); for n=1:Nx xbuffer=[x(n) xbuffer(1:end-1)]; y=[g1; g2]*xbuffer(1:Ng).'*Ts; ybuffer=[ybuffer(:,2:end) y]; subplot(312), plot(tt,ybuffer(1,:)), subplot(313), plot(tt,ybuffer(2,:)) pause(0.01), if n<Nx, cla; end end y1=xcorr(x,x1)*Ts; y1=y1([end-Nbuffer+1:end]-Ng); %correlation delayed by Nx y2=xcorr(x,x2)*Ts; y2=y2([end-Nbuffer+1:end]-Ng); subplot(312), hold on, plot(tt,y1,'m:') % only for cross-check subplot(313), hold on, plot(tt,y2,'m:') t=[0:Nx-1]*Ts; simin=[t.' x.']; % Input data for Simulink

kzbin.info/www/bejne/l17KZH6ea9-Ni5Y for Impulse Response and Transfer Function (MATLAB/Simulink Simulation) kzbin.info/www/bejne/rqS6YaakrZiloa8 for Matched Filter, Correlation, and Convolution (MATLAB/Simulink Simulation) The related files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/174315-impulse-response-and-transfer-function-with-matlab-simulink. %SS_xferfcn.m clear, clf syms s Gs=1/(s^2+2*s+2); % Transfer function (1.3.1a) w=2*pi; Us=w/(s^2+w^2); % Transformed input Ys=Gs*Us; % Transformed output y=ilaplace(Ys) % y(t)=Inverse Laplace transform of Y(s) pretty(y) t=0:0.01:10; % Time range 0~10[s] yt=eval(y); % y(t) for the time range plot(t,yt)

%sig01f06_.m % to plot several models of the unit impulse function clear, clf t=[-400:400]/200+eps; % Time range to plot the graphs on %t=[-400:400]/200; D=1; % Initialize the value of D. for i=1:5 xa=(abs(t)<D/2)/D; % Rectangular function subplot(211) plot(t,xa), axis([-2 2 -3 16.5]), hold on s1=sprintf('D=%6.4f',D); text(0.6,1/D-0.2,s1) xb=sin(pi*t/D)./(pi*t/D)/D; % Sinc function %xb=sinc(t/D)/D; % Alternatively subplot(212) plot(t,xb), axis([-2 2 -3 16.5]), hold on text(0.6,1/D-0.2,s1) shg, pause D=D/2; % Reduce D by half. end

The related m-files and slx-file can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/173945-simulink-simulation-of-an-ofdm-ieee-std-802-11a-1999. Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER The following script file be downloaded at kr.mathworks.com/matlabcentral/fileexchange/173615-matlab-simulation-of-noncoherent-fsk-passband-signaling: %sim_FSK_noncoherent_.m % simulates a noncoherent FSK passband signaling in Fig. 7.4.2 %Copyleft: Won Y. Yang, [email protected], CAU for academic use only clear, clf b=1; M=2^b; % Number of bits per symbol and modulation order ss=[0:M-1].'; Tb=1e-5; Ts=b*Tb; % Bit/Symbol durations Ns=32; T=Ts/Ns; LB=4*Ns; LBN1=LB-Ns+1; % subinterval in Tb Es=2; % the energy of signal waveform % BFSK signal waveforms wc=40*pi/Ts; dw=2*pi/Ts; w=wc+[0:M-1]*dw; wcT=wc*T; t=[0:Ns-1]*T; nd=3; % Number of delay samples td=nd*T; % Time difference between modulator and demodulator for m=1:M sw(m,:)=sqrt(2*Es/Ts)*cos(w(m)*t); %su(2*m-1:2*m,:)=sqrt(2/Ts)*[cos(w(m)*(t-td)); sin(w(m)*(t-td))]; su(2*m-1:2*m,:)=sqrt(2/Ts)*[cos(w(m)*t); sin(w(m)*t)]; end suT=su*T; SNRdBs=[1:11]; N_Iter=5000; % Range of SNRbdB and # of iterations for iter=1:length(SNRdBs) SNRbdB=SNRdBs(iter); SNRb=10^(SNRbdB/10); N0=2*(Es/b)/SNRb; sigma2=N0/2; sgmsT=sqrt(sigma2/T); sws=zeros(1,LB); r=zeros(1,LB); yr=zeros(M,LB); yrn=zeros(M,1); nose=0; % Number of symbol errors for k=1:N_Iter i=ceil(rand*M); s=ss(i,:); % Transmitted signal for n=1:Ns % Operation per symbol interval sws=[sws(2:LB) sw(i,n)]; wct=wcT*(n-1); bp_noise=randn*cos(wct)-randn*sin(wct); r=[r(2:LB) sws(end-nd)+sgmsT*bp_noise]; % Received signal possibly delayed ycss=suT*r(LBN1:LB)'; for m=1:M, yrn(m)=ycss(2*m-1)^2+ycss(2*m)^2; end yr=[yr(:,2:LB) yrn]; % Envelope detector output end %Detector(DTR) [yremax,mmax]=max(yrn); d=ss(mmax,:); nose=nose+sum(s~=d); if nose>100; break; end end pose(iter)=nose/k; end % Plot the simulated and theoretical SER curves SNRbdBt=0:0.1:SNRdBs(end); %SNRbt=10.^(SNRbdBt/10); poset=exp(-SNRbt/4)/2; %Eq.(7.2.12) poset=prob_error(SNRbdBt,'FSK',b,'SER','noncoherent'); semilogy(SNRbdBt,poset,'k', SNRdBs,pose,'b*') title(['SER for noncoherent ' num2str(M) '-FSK Signaling']) grid on, set(gca,'fontsize',9), shg %---------------------- function p=prob_error(SNRbdB,signaling,b,opt1,opt2) % Finds the symbol/bit error probability for given SNRbdB=Eb/(N0/2)[dB](Table 7.1) % Note that EbN0dB=SNRbdB-3. % Copyleft: Won Y. Yang, [email protected], CAU for academic use only if nargin<5, opt2='coherent'; end % opt2='coherent' or 'noncoherent' if nargin<4, opt1='BER'; end % opt1='SER' or 'BER' M=2^b; SNRb=10.^(SNRbdB/10); NSNR=length(SNRbdB); if signaling(1:3)=='ASK' % ASK (PAM) if lower(opt2(1))=='c' % ASK coherent p=2*(M-1)/M*Q(sqrt(3*b*SNRb/(M^2-1))); % Eq. (7.1.5) if lower(opt1(1))=='b', p = p/b; end else % ASK noncoherent if b==1, p=exp(-SNRb/4)/2; end %Eq.(7.1.22) end elseif signaling(1:3)=='FSK' tmp=M/2/(M-1); f5251_=@(x,SNRb,b)Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1); % Eq.(5.2.51) if lower(opt2(1))=='c' % FSK coherent if b==1 p=Q(sqrt(SNRb/2)); %Eq.(7.2.9) else Tol=1e-8; % tolerance used for 'quad' integration for i=1:NSNR p(i)=1-Gauss_Hermite(f5251_,10,SNRb(i),b)/sqrt(pi); end end else % FSK noncoherent for i=1:NSNR p(i)=(M-1)/2*exp(-b*SNRb(i)/4); tmp1=M-1; for m=2:M-1 tmp1=-tmp1*(M-m)/m; p(i)=p(i)+tmp1/(m+1)*exp(-m*b*SNRb(i)/2/(m+1)); % Eq.(7.2.18) end end end if lower(opt1(1))=='b'&b>1, p=p*tmp; end % Eq.(7.2.19) elseif signaling(1:3)=='PSK' p=(1+(b>1))*Q(sqrt(b*SNRb)*sin(pi/M)); %Eq.(7.3.7) if lower(opt1(1))=='b'&b>1, p=p/b; end elseif signaling(1:3)=='DPS' %DPSK if b==1 p=exp(-SNRb/2)/2; %Eq.(7.4.9) else for i=1:NSNR EbN0=SNRb(i)/2; %b; % From the MATLAB built-in function berawgn() tol=1e-6; % tol=max(min(1e-4/EbN0.^5,1e-4),eps); int_dpsk=@(y,EbN0,M,b)exp(-b*EbN0*(sin(pi/M))^2./(1+sqrt(1-(sin(pi/M)).^2).*cos(y))); %?????? p(i)=1/pi*quad(int_dpsk,1e-5,pi*(1-1/M),tol,[],EbN0,M,b); end end if lower(opt1(1))=='b', p=p/b; end elseif signaling(1:3)=='QAM' L=2^(ceil(b/2)); N = M/L; tmpL=1-2*(L-1)/L*Q(sqrt(3*b/2/(L^2-1)*SNRb)); tmpN=1-2*(N-1)/N*Q(sqrt(3*b/2/(N^2-1)*SNRb)); p=1-tmpL.*tmpN; %Eq.(7.5.6) if lower(opt1(1))=='b'&b>1, p = p/b; end end end function I=Gauss_Hermite(f,N,varargin) [t,w]=Gausshp(N); ft=feval(f,t,varargin{:}); I=w*ft'; end

great job

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER The following statements should be run before running the Simulink model where the Simulink model file can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/173775-simulink-simulation-of-coherent-4-fsk-passband-signaling: %run_before_running_FSK4_passband_coherent.m b=2; M=2^b; % Number of bits per symbol and Modulation order Ns=80; % Number of samples per symbol % This value of Ns=80 is crucial for the output signal powers of AWGN blocks, % one with a given value of SNRdB and one with the corresponding value of EbN0dB. % I don't know why. Ts=1e-5; T=Ts/Ns; % Symbol duration and Sampling interval A=sqrt(2/Ts); % Amplitude of sine/cosine wave for modulation df=1/Ts; fc=10*df; % Frequency separation and Carrier frequency dw=2*pi*df; wc=2*pi*fc; wc*T<pi, wc*T Target_no_of_error=20; Simulink_mdl='BFSK_passband_sim'; Transmitted_Signal_Power=1/Ts SNRdB=3; Received_Signal_Power_SNR=Transmitted_Signal_Power*(1+10^(-SNRdB/10)) EbN0dB=SNRdB-10*log10(b)+10*log10(Ns)-3; Received_Signal_Power_EbN0=Transmitted_Signal_Power*(1+10^(-EbN0dB/10)*Ns/b/2) Received_Signal_Power_EbN0_=Transmitted_Signal_Power*(1+10^(-(EbN0dB+3)/10)*Ns/b) % pobe_error() vs. berawgn() %SNRdB=EbN0dB+10*log10(2*b/Ns) pose=prob_error(SNRdB,'FSK',b,'sym','coherent') ber=berawgn(EbN0dB,'fsk',M,'coherent'); pose1=ber*2*(M-1)/M EbN0dB_=SNRdB-3; % This relation between SNRdB and EbN0dB fits berawgn(). ber_=berawgn(EbN0dB_,'fsk',M,'noncoherent'); pose2=ber_*2*(M-1)/M %discrepancy1=norm(pobet-pobet1)/norm(pobet)

function x=ilaplace_my(B,A) % To find the inverse Laplace transform of B(s)/A(s) using residue() % Copyleft: Won Y. Yang, [email protected], CAU for academic use only if nargin<2&numel(B)>1 for n=1:numel(B) [Num,Den]=numden(simplify(B(n))); Num=sym2poly(Num); Den=sym2poly(Den); x{n,1}=ilaplace_my(Num,Den); end return; end if ~isnumeric(B), [B,A]=numden(simplify(B)); B=sym2poly(B); A=sym2poly(A); end [r,p,k]= residue(B,A); EPS = 1e-4; N= length(r); x=[]; n=1; m=1; while n<=N %if ~isempty(x), x = [x ' + ']; end if n<N & abs(imag(p(n)))>EPS & abs(sum(imag(p([n n+1]))))<EPS sigma=real(p(n)); w=imag(p(n)); Kc=2*real(r(n)); Ks=-2*imag(r(n)); sigma_=num2str(sigma); w_=num2str(w); Kc_=num2str(Kc); Ks_=num2str(Ks); if Kc>0&~isempty(x), x = [x ' + ']; end if abs(sigma)>EPS, x = [x 'exp(' sigma_ '*t).*']; end if abs(sigma)>EPS&abs(Kc)>EPS&abs(Ks)>EPS, x=[x '(']; end if abs(Kc)>EPS if Kc~=1, x=[x Kc_ '*']; end x=[x 'cos(' w_ '*t)']; end if abs(Ks)>EPS if Ks>EPS&abs(Kc)>EPS, x = [x '+']; end if abs(Ks-1)>EPS if Ks==-1, x=[x '-']; else x=[x Ks_ '*']; end end x=[x 'sin(' w_ '*t)']; if abs(sigma)>EPS&abs(Kc)>EPS, x=[x ') ']; end end n = n + 2; elseif n<=N & abs(imag(r(n)))<EPS if n>1, if p(n)==p(n-1), m=m+1; else m=1; end; end if abs(r(n))>EPS if r(n)>0&~isempty(x), x = [x '+ ']; end r(n)=r(n)/factorial(m-1); if abs(r(n)-1)>EPS, x=[x num2str(r(n)) '*']; end if abs(r(n)+1)<EPS, x=[x ' - ']; end if m==2, x = [x 't.*']; elseif m>2, x = [x 't.^' num2str(m-1) '.*']; end if abs(p(n))>EPS, x=[x 'exp(' num2str(p(n)) '*t) ']; else x=[x '1']; end end n = n+1; else fprintf('ilaplace_my() cannot deal with X(s) having multiple complex poles '); break; end end if ~isempty(k) Kd = num2str(k(end)); if k(end)>0, x = [x ' + ']; end x = [x Kd '*dirac(t)']; end %Usage 1: B=[4 0 2], A=[6 4 4 1]; ilaplace_my(B,A) %Usage 2: syms s; ilaplace_my((4*s^2+2)/(6*s^3+4*s^2+4*s+1))

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER %sim_DPSK_passband_.m % simulates a digital communication system with QDPSK signal waveforms in Fig. 7.9 clear, clf b=2; M=2^b; % Number of bits per symbol and Modulation order ss=bin2gray(0:M-1,'psk',M); % Set of Gray-coded b-bit symbols %ss=[0 0; 0 1; 1 1; 1 0]; %[gc_int,ss]=gray_codes(b); Tb=1e-5; Ts=b*Tb; % Bit time & Symbol time Nb=32; Ns=b*Nb; % Number of subintervals in Ts T=Ts/Ns; LB=4*Ns; LBN1=LB-Ns+1; % Sample time and Buffer size Es=2; sqEs=sqrt(Es); % Energy of signal waveform % Carrier frequency wc=8*pi/Ts; t=[0:Ns-1]*T; wcT=wc*T; phases=[0:M-1]*2*pi/M; gs='>^<vs*x+d'; % PSK signal waveforms for m=1:M sw(m,:)=sqrt(2*Es/Ts)*cos(wc*t+phases(m)); % Signal waveforms to transmit end su=sqrt(2/Ts)*[cos(wc*t); -sin(wc*t)]; suT=su*T; nd=1; % Number of delay samples SNRdBs=[1:11]; N_Iter=5e5; % Range of SNRbdB and Number of iterations for iter=1:length(SNRdBs) SNRbdB=SNRdBs(iter); SNRb=10^(SNRbdB/10); sigma2=(Es/b)/SNRb; sgmsT=sqrt(sigma2/T); sws=zeros(1,LB); yr=zeros(2,LB); nobe=0; %number of bit errors to see how well DPSK combats the delay is0=0; % Initial signal th0=phases(3); % Initial guess (possibly wrong) Discrepancy=0; for k=1:N_Iter i=floor(rand*M); s=ss(i+1,:); %transmitted information is=mod(is0+i,M); is1=is+1; %transmitted signal for n=1:Ns %operation per symbol interval sws=[sws(2:LB) sw(is1,n)]; wct=wcT*(n-1); bp_noise=randn*cos(wct)-randn*sin(wct); rn=sws(end-nd)+sgmsT*bp_noise; % A sample of received signal yr=[yr(:,2:LB) suT(:,n)*rn]; % Multiplication for correlator end ycsk=sum(yr(:,LBN1:LB)'); % Sampled quadrature correlator outputs - DTR input %Detector(DTR) th=atan2(ycsk(2),ycsk(1)); dth=th-th0; [dmin,mmin]=min(abs(exp(j*dth)-exp(j*phases))); d=ss(mmin,:); nobe=nobe+sum(s~=d); if nobe>100, break; end is0=is; th0=th; %update the previous signal and theta end %Discrepancy pobe(iter)=nobe/(k*b); end % Plot the simulated and theoretical SER curves. SNRbdBt=0:0.1:SNRdBs(end); EbN0dBt=SNRbdBt-3; [BERt,SERt]=berawgn(EbN0dBt,'dpsk',M); semilogy(SNRbdBt,BERt,'b', SNRdBs,pobe,'r*') title(['BER for passband ' num2str(M) '-DPSK Signaling']) set(gca,'fontsize',9), shg %----------------------- function [grayCode,grayCode_s,map]=bin2gray(xin,signaling,M) % Input: % xin = Integer(s) to be graycoded. % signaling = 'ask'|'psk'|'qam' % M = Modulation order % Output: % grayCode = Binary form % grayCode_s = String form % map = Grade code map in string form signaling=lower(signaling); if signaling=='ask'|signaling=='pam' log2M=log2(M); binVal = de2bi(xin,log2M,'left-msb'); % Convert binary to Gray code grayCode = binVal; grayCode(:,2:end) = xor(binVal(:,1:end-1), binVal(:, 2:end)); binVals = de2bi([0:M-1],log2M,'left-msb'); grayCodes = binVals; grayCodes(:,2:end) = xor(binVals(:,1:end-1), binVals(:, 2:end)); for m=1:size(grayCode,1) grayCode_s(m,:)=removeBlanks(num2str(grayCode(m,:))); % String form end for m=1:size(grayCodes,1) map{m}=removeBlanks(num2str(grayCodes(m,:))); % String form end for m=1:numel(xin) grayCode_s1(m,:)=map{xin(m)+1}; % String form end elseif signaling=='qam' b=log2(M); L=2^(b/2); graycode=bin2gray([0:L-1],'ask',L); i=0; for m=1:L for n=1:L i = i+1; grayCodes(i,:)=[graycode(n,:) graycode(m,:)]; map{L-m+1,n}=removeBlanks(num2str(grayCodes(i,:))); end end %grayCode0 grayCode=grayCodes(xin+1,:); % Binary integer form for m=1:size(grayCode,1) grayCode_s(m,:)=removeBlanks(num2str(grayCode(m,:))); % String form end elseif signaling=='psk' b=ceil(log2(M)); grayCodes=[0; 1]; for n=2:b m=size(grayCodes,1); grayCodes=[zeros(m,1) grayCodes; ones(m,1) flipud(grayCodes)]; end grayCode=grayCodes(xin+1,:); % Binary integer form for m=1:size(grayCodes,1) map{m}=removeBlanks(num2str(grayCodes(m,:))); % String form end for m=1:numel(xin) grayCode_s(m,:)=map{xin(m)+1}; % String form end else fprintf(" Please specify the 2nd input argument signaling as one of 'ask'|'psk','qam'") end end

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER The Simulink model file is available at kr.mathworks.com/matlabcentral/fileexchange/173550-simulink-simulation-of-dpsk-passband-signaling. %do_DPSK_passband_sim.m clear, clf b=2; M=2^b; % Number of bits per symbol and Modulation order Ns=160; % This is crucial to getting a realistic BER performance Ts=1e-5; T=Ts/Ns; % Symbol duration and Sampling interval A=sqrt(2/Ts); wc=2*pi*10/Ts; % Carrier frequency[rad/s] d=5; % Some delay to see how well the DPSK combats a delay Target_no_of_error=100; Transmitted_Signal_Power=1/Ts; Simulink_mdl='DPSK_passband_sim'; SNRdBs=[1:2:11]; Received_Signal_Power=Transmitted_Signal_Power*(1+10^(-SNRdBs(end)/10)*Ts/b/T) for i=1:length(SNRdBs) SNRdB=SNRdBs(i); EbN0dB=SNRdB-3; %10*log10(b); out=sim(Simulink_mdl,2e5*Ts); SERs(i)=ser(end,1); fprintf('%d-ary %s with SNRdB=%5.2f[dB] yields SER=%10.4e ',M,Simulink_mdl,SNRdB,SERs(i)) end SERs SNRdBt=[0:0.1:SNRdBs(end)]; EbN0dBt=SNRdBt-3; poset=prob_error(SNRdBt,'DPSK',b,'sym'); [BERt,SERt]=berawgn(EbN0dBt,'DPSK',M); %SERt=b*BERt; semilogy(SNRdBt,poset,'k', SNRdBs,SERs,'*', SNRdBt,SERt,'m:') title(sprintf('SER for %2d-DPSK passband signaling',M)) xlabel('SNR[dB]'), shg function mo=detector_DPSK(dth,M) phases=2*pi/M*[0:M-1]; [dmin,mmin]=min(abs(exp(j*dth)-exp(j*phases))); mo=mmin-1; function p=prob_error(SNRbdB,signaling,b,opt1,opt2) % Finds the symbol/bit error probability for given SNRbdB=Eb/(N0/2)[dB](Table 7.1) % Note that EbN0dB=SNRbdB-3. % Copyleft: Won Y. Yang, [email protected], CAU for academic use only if nargin<5, opt2='coherent'; end % opt2='coherent' or 'noncoherent' if nargin<4, opt1='BER'; end % opt1='SER' or 'BER' M=2^b; SNRb=10.^(SNRbdB/10); NSNR=length(SNRbdB); if signaling(1:3)=='ASK' % ASK (PAM) if lower(opt2(1))=='c' % ASK coherent p=2*(M-1)/M*Q(sqrt(3*b*SNRb/(M^2-1))); % Eq. (7.1.5) if lower(opt1(1))=='b', p = p/b; end else % ASK noncoherent if b==1, p=exp(-SNRb/4)/2; end %Eq.(7.1.22) end elseif signaling(1:3)=='FSK' tmp=M/2/(M-1); f5251_=@(x,SNRb,b)Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1); % Eq.(5.2.51) if lower(opt2(1))=='c' % FSK coherent if b==1 p=Q(sqrt(SNRb/2)); %Eq.(7.2.9) else Tol=1e-8; % tolerance used for 'quad' integration for i=1:NSNR p(i)=1-Gauss_Hermite(f5251_,10,SNRb(i),b)/sqrt(pi); end end else % FSK noncoherent for i=1:NSNR p(i)=(M-1)/2*exp(-b*SNRb(i)/4); tmp1=M-1; for m=2:M-1 tmp1=-tmp1*(M-m)/m; p(i)=p(i)+tmp1/(m+1)*exp(-m*b*SNRb(i)/2/(m+1)); % Eq.(7.2.18) end end end if lower(opt1(1))=='b'&b>1, p=p*tmp; end % Eq.(7.2.19) elseif signaling(1:3)=='PSK' p=(1+(b>1))*Q(sqrt(b*SNRb)*sin(pi/M)); %Eq.(7.3.7) if lower(opt1(1))=='b'&b>1, p=p/b; end elseif signaling(1:3)=='DPS' %DPSK if b==1 p=exp(-SNRb/2)/2; %Eq.(7.4.9) else for i=1:NSNR EbN0=SNRb(i)/2; %b; % From the MATLAB built-in function berawgn() tol=1e-6; % tol=max(min(1e-4/EbN0.^5,1e-4),eps); int_dpsk=@(y,EbN0,M,b)exp(-b*EbN0*(sin(pi/M))^2./(1+sqrt(1-(sin(pi/M)).^2).*cos(y))); %?????? p(i)=1/pi*quad(int_dpsk,1e-5,pi*(1-1/M),tol,[],EbN0,M,b); end end if lower(opt1(1))=='b', p=p/b; end elseif signaling(1:3)=='QAM' L = 2^(ceil(b/2)); N = M/L; tmpL = 1-2*(L-1)/L*Q(sqrt(3*b/2/(L^2-1)*SNRb)); tmpN = 1-2*(N-1)/N*Q(sqrt(3*b/2/(N^2-1)*SNRb)); p = 1-tmpL.*tmpN; %Eq.(7.5.6) if lower(opt1(1))=='b'&b>1, p = p/b; end end end function I=Gauss_Hermite(f,N,varargin) [t,w]=Gausshp(N); ft=feval(f,t,varargin{:}); I=w*ft'; end

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER The Simulink model file is available at kr.mathworks.com/matlabcentral/fileexchange/173040-simulation-of-qam-passband-signaling-with-simulink.

Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulation of noncoherent FSK passband signaling to measure the SER with Simulink. Watch kzbin.info/www/bejne/rWW0iZhoeMyXfKc for Simulation of noncoherent FSK passband signaling to measure the SER with MATLAB. Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for Simulation of coherent FSK passband signaling to measure the SER with MATLAB. Watch kzbin.info/www/bejne/apS0dJufadNki7s for Simulation of QAM passband signaling to measure the BER with MATLAB Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulation of QAM passband signaling to measure the SER with Simulink Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for Simulation of PSK passband signaling to measure the BER with MATLAB Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulation of PSK passband signaling to measure the BER with Simulink

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER The MATLAB script file can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/173590-matlab-simulation-of-psk-passband-signaling-to-measure-ber. %sim_PSK_passband.m % simulates a digital communication system with QPSK signal waveforms in Fig.7.7 %Copyleft: Won Y. Yang, [email protected], CAU for academic use only clear, clf b=2; M=2^b; % Number of bits per symbol and Modulation order ss=bin2gray([0:M-1],'psk',M); % Symbol (Alphabet) set Tb=1e-5; Ts=b*Tb; % Bit/Symbol durations Nb=32; Ns=b*Nb; % Number of samples per Ts and Sample time T=Ts/Ns; LB=4*Ns; LBN1=LB-Ns+1; % Sample time, Buffer size Es=2; sqEs=sqrt(Es); % Energy of signal waveform wc=8*pi/Ts; % Carrier frequency t=[0:Ns-1]*T; wcT=wc*T; nd=0; % Number of delay samples % PSK signal waveforms phases=2*pi/M*[0:M-1]; gs='>^<vs*x+do>^<vs*x+d'; for m=1:M sw(m,:)=sqrt(2*Es/Ts)*cos(wc*t+phases(m)); end su=sqrt(2/Ts)*[cos(wc*t); -sin(wc*t)]; suT=su*T; SNRbdBs=[1:10]; % Range of SNRbdB N_Iter=10000; % Number of iterations for getting the probability of error for iter=1:length(SNRbdBs) SNRbdB= SNRbdBs(iter); SNRb=10^(SNRbdB/10); N0=2*(Es/b)/SNRb; sigma2=N0/2; sgmsT=sqrt(sigma2/T); sws=zeros(1,LB); r=zeros(1,LB); yr=zeros(2,LB); nobe=0; %number of bit errors for k=1:N_Iter % Loop iterations for each SNRdB value i=ceil(rand*M); s=ss(i,:); % Transmitted signal for n=1:Ns % Operation per symbol interval sws=[sws(2:LB) sw(i,n)]; %Received signal wct=wcT*(n-1); bp_noise=randn*cos(wct)-randn*sin(wct); rn=sws(end-nd)+sgmsT*bp_noise; % A sample of received signal yr=[yr(:,2:LB) suT(:,n)*rn]; % Multiplication for correlator end ycsk=sum(yr(:,LBN1:LB)')'; % Sampled quadrature correlator outputs - DTR input %Detector(DTR) if M<17&iter==9&k<300 % Signal constellation diagram subplot(221), hold on plot(ycsk(1),ycsk(2),['b' gs(i)],'Markersize',5) end th=atan2(ycsk(2),ycsk(1)); if th<-pi/M, th=th+2*pi; end [thdmin,mo]=min(abs(th-phases)); d=ss(mo,:); nobe=nobe+sum(s~=d); if nobe>100, break; end end pobe(iter)=nobe/(k*b); end plot(sqEs*exp(j*phases),'ro','Markersize',6) axis([-2.5 2.5 -2.5 2.5]), axis('equal') set(gca,'fontsize',9) subplot(222) SNRbdBt=0:0.1:SNRbdBs(end); SNRbt=10.^(SNRbdBt/10); pobet=prob_error(SNRbdBt,'PSK',b); semilogy(SNRbdBt,pobet,'k', SNRbdBs,pobe,'b*') title(sprintf('BER for %2d-PSK Signaling',M)) grid on, shg %----------------------- function [grayCode,grayCode_s,map]=bin2gray(xin,signaling,M) % Input: % xin = Integer(s) to be graycoded. % signaling = 'ask'|'psk'|'qam' % M = Modulation order % Output: % grayCode = Binary form % grayCode_s = String form % map = Grade code map in string form signaling=lower(signaling); if signaling=='ask'|signaling=='pam' log2M=log2(M); binVal=de2bi(xin,log2M,'left-msb'); % Convert binary to Gray code grayCode=binVal; grayCode(:,2:end)=xor(binVal(:,1:end-1), binVal(:, 2:end)); binVals=de2bi([0:M-1],log2M,'left-msb'); grayCodes=binVals; grayCodes(:,2:end)=xor(binVals(:,1:end-1), binVals(:, 2:end)); for m=1:size(grayCode,1) grayCode_s(m,:)=removeBlanks(num2str(grayCode(m,:))); % String form end for m=1:size(grayCodes,1) map{m}=removeBlanks(num2str(grayCodes(m,:))); % String form end for m=1:numel(xin) grayCode_s1(m,:)=map{xin(m)+1}; % String form end elseif signaling=='qam' b=log2(M); L=2^(b/2); graycode=bin2gray([0:L-1],'ask',L); i=0; for m=1:L for n=1:L i = i+1; grayCodes(i,:)=[graycode(n,:) graycode(m,:)]; map{L-m+1,n}=removeBlanks(num2str(grayCodes(i,:))); end end %grayCode0 grayCode=grayCodes(xin+1,:); % Binary integer form for m=1:size(grayCode,1) grayCode_s(m,:)=removeBlanks(num2str(grayCode(m,:))); % String form end elseif signaling=='psk' b=ceil(log2(M)); grayCodes=[0; 1]; for n=2:b m=size(grayCodes,1); grayCodes=[zeros(m,1) grayCodes; ones(m,1) flipud(grayCodes)]; end grayCode=grayCodes(xin+1,:); % Binary integer form for m=1:size(grayCodes,1) map{m}=removeBlanks(num2str(grayCodes(m,:))); % String form end for m=1:numel(xin) grayCode_s(m,:)=map{xin(m)+1}; % String form end else fprintf(" Please specify the 2nd input argument signaling as one of 'ask'|'psk','qam'") end end function p=prob_error(SNRbdB,signaling,b,opt1,opt2) % Finds the symbol/bit error probability for given SNRbdB=Eb/(N0/2)[dB](Table 7.1) % Note that EbN0dB=SNRbdB-3. % Copyleft: Won Y. Yang, [email protected], CAU for academic use only if nargin<5, opt2='coherent'; end % opt2='coherent' or 'noncoherent' if nargin<4, opt1='BER'; end % opt1='SER' or 'BER' M=2^b; SNRb=10.^(SNRbdB/10); NSNR=length(SNRbdB); if signaling(1:3)=='ASK' % ASK (PAM) if lower(opt2(1))=='c' % ASK coherent p=2*(M-1)/M*Q(sqrt(3*b*SNRb/(M^2-1))); % Eq. (7.1.5) if lower(opt1(1))=='b', p = p/b; end else % ASK noncoherent if b==1, p=exp(-SNRb/4)/2; end %Eq.(7.1.22) end elseif signaling(1:3)=='FSK' tmp=M/2/(M-1); f5251_=@(x,SNRb,b)Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1); % Eq.(5.2.51) %f7208=inline('Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1).*exp(-x.^2)','x','SNRb','b'); if lower(opt2(1))=='c' % FSK coherent if b==1 p=Q(sqrt(SNRb/2)); %Eq.(7.2.9) else Tol=1e-8; % tolerance used for 'quad' integration for i=1:NSNR p(i)=1-Gauss_Hermite(f5251_,10,SNRb(i),b)/sqrt(pi); %p(i)=1-adapt_smpsn('f7208',-20,20,Tol,SNRb(i),b)/sqrt(pi); % Eq.(7.2.8) %p(i)=1-quadl('f7208',-20,20,Tol,[],SNRb(i),b)/sqrt(pi); end end else % FSK noncoherent for i=1:NSNR p(i)=(M-1)/2*exp(-b*SNRb(i)/4); tmp1=M-1; for m=2:M-1 tmp1=-tmp1*(M-m)/m; p(i)=p(i)+tmp1/(m+1)*exp(-m*b*SNRb(i)/2/(m+1)); % Eq.(7.2.18) end end end if lower(opt1(1))=='b'&b>1, p=p*tmp; end % Eq.(7.2.19) elseif signaling(1:3)=='PSK' p=(1+(b>1))*Q(sqrt(b*SNRb)*sin(pi/M)); %Eq.(7.3.7) if lower(opt1(1))=='b'&b>1, p=p/b; end elseif signaling(1:3)=='DPS' %DPSK if b==1 p=exp(-SNRb/2)/2; %Eq.(7.4.9) else for i=1:NSNR EbN0=SNRb(i)/2; %b; % From the MATLAB built-in function berawgn() tol=1e-6; % tol=max(min(1e-4/EbN0.^5,1e-4),eps); int_dpsk=@(y,EbN0,M,b)exp(-b*EbN0*(sin(pi/M))^2./(1+sqrt(1-(sin(pi/M)).^2).*cos(y))); %?????? p(i)=1/pi*quad(int_dpsk,1e-5,pi*(1-1/M),tol,[],EbN0,M,b); end end if lower(opt1(1))=='b', p=p/b; end elseif signaling(1:3)=='QAM' L=2^(ceil(b/2)); N = M/L; tmpL=1-2*(L-1)/L*Q(sqrt(3*b/2/(L^2-1)*SNRb)); tmpN=1-2*(N-1)/N*Q(sqrt(3*b/2/(N^2-1)*SNRb)); p=1-tmpL.*tmpN; %Eq.(7.5.6) if lower(opt1(1))=='b'&b>1, p = p/b; end end end function I=Gauss_Hermite(f,N,varargin) [t,w]=Gausshp(N); ft=feval(f,t,varargin{:}); I=w*ft'; end

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER %sim_FSK_coherent.m % simulates a coherent FSK passband signaling in Fig.7.4.1 %Copyleft: Won Y. Yang, [email protected], CAU for academic use only clear, clf KC=2; % 0/1/2 no time offset/time offset and no sync/time offset and sync b=1; M=2^b; % Number of bits per symbol and Modulation order ss=[0:M-1].'; % Symbol (Alphabet) set Tb=1e-5; Ts=b*Tb; % Bit/Symbol durations Ns=32; T=Ts/Ns; % Number of samples per Ts LB=4*Ns; LBN1=LB-Ns+1; % Number of samples kept in buffer corresponding 4 symbol durations Es=2; % the energy of signal waveform % Carrier frequency dw=2*pi/Ts; wc=8*dw; w=wc+[0:M-1]*dw; % M tone frequencies wcT=wc*T; t=[0:Ns-1]*T; % Time difference & Synchronization between modulator/demodulator if KC<1, nd=0; else nd=1; end % Number of delay samples if KC<2, td=0; else td=nd*T; end % for synchronization % FSK signal waveforms and Correlator signal waveforms for detection with some time difference for m=1:M sw(m,:)=sqrt(2*Es/Ts)*cos(w(m)*t); % Transmit signal waveforms su(m,:)=sqrt(2/Ts)*cos(w(m)*(t-td)); % Eq.(7.2.11) end suT=su*T; % Signal waveforms used in the correlator for detection SNRdBs=[1:11]; N_Iter=1e5; % Range of SNRbdB and # of iterations for iter=1:length(SNRdBs) SNRbdB=SNRdBs(iter); SNRb=10^(SNRbdB/10); N0=2*(Es/b)/SNRb; sigma2=N0/2; sgmsT=sqrt(sigma2/T); sws=zeros(1,LB); r=zeros(1,LB); yr=zeros(M,LB); nose=0; % Number of symbol errors for k=1:N_Iter i=ceil(rand*M); s=ss(i,:); % Transmitted signal for n=1:Ns % Operation per symbol interval sws=[sws(2:LB) sw(i,n)]; % i-th wave form (i: random integer among [0:M-1] wct=wcT*(n-1); bp_noise=randn*cos(wct)-randn*sin(wct); r=[r(2:LB) sws(end-nd)+sgmsT*bp_noise]; %received signal possibly delayed yrn=suT*r(LBN1:LB).'; yr=[yr(:,2:LB) yrn]; % Correlator output - DTR input end %Detector(DTR) [yremax,mmax]=max(yrn); d=ss(mmax,:); nose=nose+sum(s~=d); if nose>100; break; end end pose(iter)=nose/k; end SNRbdBt=0:0.1:SNRdBs(end); poset=prob_error(SNRbdBt,'FSK',b,'SER'); semilogy(SNRbdBt,poset,'k', SNRdBs,pose,'b*') ylim([0.005 1]), grid on title(['SER for coherent ' num2str(M) '-FSK Signaling']) grid on, shg %----------------------- function p=prob_error(SNRbdB,signaling,b,opt1,opt2) % Finds the symbol/bit error probability for given SNRbdB=Eb/(N0/2)[dB](Table 7.1) % Note that EbN0dB=SNRbdB-3. % Copyleft: Won Y. Yang, [email protected], CAU for academic use only if nargin<5, opt2='coherent'; end % opt2='coherent' or 'noncoherent' if nargin<4, opt1='BER'; end % opt1='SER' or 'BER' M=2^b; SNRb=10.^(SNRbdB/10); NSNR=length(SNRbdB); if signaling(1:3)=='ASK' % ASK (PAM) if lower(opt2(1))=='c' % ASK coherent p=2*(M-1)/M*Q(sqrt(3*b*SNRb/(M^2-1))); % Eq. (7.1.5) if lower(opt1(1))=='b', p = p/b; end else % ASK noncoherent if b==1, p=exp(-SNRb/4)/2; end %Eq.(7.1.22) end elseif signaling(1:3)=='FSK' tmp=M/2/(M-1); f5251_=@(x,SNRb,b)Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1); % Eq.(5.2.51) %f7208=inline('Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1).*exp(-x.^2)','x','SNRb','b'); if lower(opt2(1))=='c' % FSK coherent if b==1 p=Q(sqrt(SNRb/2)); %Eq.(7.2.9) else Tol=1e-8; % tolerance used for 'quad' integration for i=1:NSNR p(i)=1-Gauss_Hermite(f5251_,10,SNRb(i),b)/sqrt(pi); %p(i)=1-adapt_smpsn('f7208',-20,20,Tol,SNRb(i),b)/sqrt(pi); % Eq.(7.2.8) %p(i)=1-quadl('f7208',-20,20,Tol,[],SNRb(i),b)/sqrt(pi); end end else % FSK noncoherent for i=1:NSNR p(i)=(M-1)/2*exp(-b*SNRb(i)/4); tmp1=M-1; for m=2:M-1 tmp1=-tmp1*(M-m)/m; p(i)=p(i)+tmp1/(m+1)*exp(-m*b*SNRb(i)/2/(m+1)); % Eq.(7.2.18) end end end if lower(opt1(1))=='b'&b>1, p=p*tmp; end % Eq.(7.2.19) elseif signaling(1:3)=='PSK' p=(1+(b>1))*Q(sqrt(b*SNRb)*sin(pi/M)); %Eq.(7.3.7) if lower(opt1(1))=='b'&b>1, p=p/b; end elseif signaling(1:3)=='DPS' %DPSK for i=1:NSNR EbN0=SNRb(i)/b; % From the MATLAB built-in function berawgn() tol=1e-6; % tol=max(min(1e-4/EbN0.^5,1e-4),eps); int_dpsk=@(y,EbN0,M,b)exp(-b*EbN0*(sin(pi/M))^2./(1+sqrt(1-(sin(pi/M)).^2).*cos(y))); %?????? p(i)=1/pi*quad(int_dpsk,1e-5,pi*(1-1/M),tol,[],EbN0,M,b); end if lower(opt1(1))=='b', p=p/b; end elseif signaling(1:3)=='QAM' L = 2^(ceil(b/2)); N = M/L; tmpL = 1-2*(L-1)/L*Q(sqrt(3*b/2/(L^2-1)*SNRb)); tmpN = 1-2*(N-1)/N*Q(sqrt(3*b/2/(N^2-1)*SNRb)); p = 1-tmpL.*tmpN; %Eq.(7.5.6) if lower(opt1(1))=='b'&b>1, p = p/b; end end end function I=Gauss_Hermite(f,N,varargin) [t,w]=Gausshp(N); ft=feval(f,t,varargin{:}); I=w*ft'; end

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER The following statements should be run before running the Simulink model where the Simulink model file can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/173715-simulink-simulation-of-noncoherent-bfsk-passband-signaling: %run_before_running_BFSK_passband_noncoherent.m b=1; M=2^b; % Number of bits per symbol and Modulation order Ns=80; % Number of samples per symbol % This value of Ns=80 is crucial for the output signal powers of AWGN blocks, % one with a given value of SNRdB and one with the corresponding value of EbN0dB. % I don't know why. Ts=1e-5; T=Ts/Ns; % Symbol duration and Sampling interval A=sqrt(2/Ts); % Amplitude of sine/cosine wave for modulation df=1/Ts; fc=10*df; % Frequency separation and Carrier frequency dw=2*pi*df; wc=2*pi*fc; wc*T<pi, wc*T Target_no_of_error=20; Simulink_mdl='BFSK_passband_sim'; Transmitted_Signal_Power=1/Ts SNRdB=3; Received_Signal_Power_SNR=Transmitted_Signal_Power*(1+10^(-SNRdB/10)) EbN0dB=SNRdB-10*log10(b)+10*log10(Ns)-3; Received_Signal_Power_EbN0=Transmitted_Signal_Power*(1+10^(-EbN0dB/10)*Ns/b/2) Received_Signal_Power_EbN0_=Transmitted_Signal_Power*(1+10^(-(EbN0dB+3)/10)*Ns/b)

Watch kzbin.info/www/bejne/hXiqo4FnitqXhMU for Simulink Simulation of OFDM communication system to get its BER Watch kzbin.info/www/bejne/Z4vEXphpd8qhn9U for Simulink Simulation of coherent 4-FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/a568gKakaMx9i9k for MATLAB Simulation of DPSK (Differential PSK) passband signaling to measure the BER Watch kzbin.info/www/bejne/eIeseZyhataioa8 for Simulink Simulation of DPSK (Differential PSK) passband signaling to measure the SER Watch kzbin.info/www/bejne/aoW3eoufjbiIb8U for Simulink Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/h3a9nHqwg9yWfMk for MATLAB Simulation of noncoherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/bXXbfGlpZrCqq80 for MATLAB Simulation of coherent FSK passband signaling to measure the SER Watch kzbin.info/www/bejne/apS0dJufadNki7s for MATLAB Simulation of QAM passband signaling to measure the BER Watch kzbin.info/www/bejne/gZWvoICfYst8i5Y for Simulink Simulation of QAM passband signaling to measure the SER Watch kzbin.info/www/bejne/nJu3nJWEnd51nZo for MATLAB Simulation of PSK passband signaling to measure the BER Watch kzbin.info/www/bejne/h5u0hZJmqtFqb80 for Simulink Simulation of PSK passband signaling to measure the BER %sim_QAM_passband.m % simulates a digital communication system in Fig.7.13 % with QAM signal waveforms in Fig.7.11 clear, clf b=4; M=2^b; L=2^(b/2); % # of bits per symbol and modulation order Tb=1e-5; Ts=b*Tb; % Bit/Symbol duration Nb=16; Ns=b*Nb; % # of subintervals in Ts T=Ts/Ns; LB=4*Ns; LBN1=LB-Ns+1; % the subinterval length and buffer size ssc=bin2gray(0:L-1,'ask',L); sss=ssc; wc=8*pi/Ts; wcT=wc*T; t=[0:Ns-1]*T; % Unit signal waveforms for the quadrature multiplier of transmitter su=sqrt(2/Ts)*[cos(wc*t); -sin(wc*t)]; suT=su*T; % Unit signal waveforms for the quadrature multiplier of receiver %Esum=0; % 16-QAM signal waveforms corresponding to rectangular constellation for i=1:L for l=1:L s(i,l,1)=2*i-L-1; s(i,l,2)=2*l-L-1; %In-phase/quadrature amplitude %Esum= Esum+s(i,l,1)^2 +s(i,l,2)^2; ss(L*(l-1)+i,:)=[ssc(i,:) sss(l,:)]; sw(L*(l-1)+i,:)=s(i,l,1)*su(1,:)+s(i,l,2)*su(2,:); end end %Eav_=Esum/M Eav=2*(M-1)/3 % Eq.(7.5.4a) Es=2; % the energy of signal waveform A=sqrt(Es/Eav); sw=A*sw; levels=A*[-(L-1):2:L-1]; gs=['s+xs'; 'd<^d'; '*v>*'; 's+xs']; % graphic symbols SNRdBs=[1:15]; N_Iter=10000; % Range of SNRbdB and # of iterations for iter=1:length(SNRdBs) SNRbdB=SNRdBs(iter); SNRb=10^(SNRbdB/10); %SNRdB=SNRbdB+10*log10(b); sigma2=(Es/b)/SNRb; sgmsT=sqrt(sigma2/T); yr=zeros(2,LB); % Initialize the multiplier buffer nobe=0; % Initialize the number of bit errors to be accumulated for k=1:N_Iter im=ceil(rand*L); in=ceil(rand*L); imn=(in-1)*L+im; s=ss(imn,:); % Transmitted signal for n=1:Ns % Operation per symbol interval wct=wcT*(n-1); bp_noise=randn*cos(wct)-randn*sin(wct); % Bandpass noise rn=sw(imn,n)+sgmsT*bp_noise; % Received signal with noise %rn= awgn(sw(imn,n),SNRdB-10*log10(Ns/Ts)-3); %received signal yr=[yr(:,2:LB) suT(:,n)*rn]; % Multiplier end ycsk=sum(yr(:,LBN1:LB)'); % Sampled correlator output - DTR input %yck=ycsk(1); ysk=ycsk(2); %Detector(DTR)levels if iter==10&k<300 %signal constellation diagram subplot(221), hold on plot(ycsk(1),ycsk(2),['b' gs(L+1-in,im)],'Markersize',5) end % Indices of the in-phase/quadrature levels closest to the correlator outputs [dmin_i,mi]=min(abs(ycsk(1)-levels)); [dmin_l,ml]=min(abs(ycsk(2)-levels)); %fprintf(' yre=%7.2f %7.2f, im=%2d, in=%2d, mi=%2d, ml=%2d', yre(1),yre(2),im,in,mi,ml) d=ss((ml-1)*L+mi,:); % Detected signal nobe=nobe+sum(s~=d); % Number of bit errors if nobe>100; break; end end pobe(iter)= nobe/(k*b); end % Add the reference signal points to the signal constellation diagram for n=1:L, plot(levels,levels(n)*ones(1,L),'ro'); end axis([-2 2 -2 2]), axis('equal'), set(gca,'fontsize',9) title(['Received signal constellation for SNRdB=' num2str(SNRdBs(10))]) % Plot the BERs obtained from simulation with the theoretical BER curve SNRbdBt=0:.1:15; SNRbt=10.^(SNRbdBt/10); %Pm=2*(1-1/L)*Q(sqrt(3/2*b*SNRbt/(M-1))); % Eq.(7.1.5) %poe_on_theory= (1-(1-Pm).^2)/b; % Eq.(7.5.8) poe_on_theory=prob_error(SNRbdBt,'QAM',b,'BER'); subplot(222) semilogy(SNRbdBt,poe_on_theory,'k-', SNRdBs,pobe,'b*') axis([0 SNRdBs(end) 1e-4 1]), set(gca,'fontsize',9), grid on title('Probability of Bit Error for (16-ary) QAM Signaling') set(gcf,'color','white') function [grayCode,grayCode_s,map]=bin2gray(xin,signaling,M) % Input: % xin = Integer(s) to be graycoded. % signaling = 'ask'|'psk'|'qam' % M = Modulation order % Output: % grayCode = Binary form % grayCode_s = String form % map = Grade code map in string form signaling=lower(signaling); if signaling=='ask'|signaling=='pam' log2M=log2(M); binVal = de2bi(xin,log2M,'left-msb'); % Convert binary to Gray code grayCode = binVal; grayCode(:,2:end) = xor(binVal(:,1:end-1), binVal(:, 2:end)); binVals = de2bi([0:M-1],log2M,'left-msb'); grayCodes = binVals; grayCodes(:,2:end) = xor(binVals(:,1:end-1), binVals(:, 2:end)); for m=1:size(grayCode,1) grayCode_s(m,:)=removeBlanks(num2str(grayCode(m,:))); % String form end for m=1:size(grayCodes,1) map{m}=removeBlanks(num2str(grayCodes(m,:))); % String form end for m=1:numel(xin) grayCode_s1(m,:)=map{xin(m)+1}; % String form end elseif signaling=='qam' b=log2(M); L=2^(b/2); graycode=bin2gray([0:L-1],'ask',L); i=0; for m=1:L for n=1:L i = i+1; grayCodes(i,:)=[graycode(n,:) graycode(m,:)]; map{L-m+1,n}=removeBlanks(num2str(grayCodes(i,:))); end end %grayCode0 grayCode=grayCodes(xin+1,:); % Binary integer form for m=1:size(grayCode,1) grayCode_s(m,:)=removeBlanks(num2str(grayCode(m,:))); % String form end elseif signaling=='psk' b=ceil(log2(M)); grayCodes=[0; 1]; for n=2:b m=size(grayCodes,1); grayCodes=[zeros(m,1) grayCodes; ones(m,1) flipud(grayCodes)]; end grayCode=grayCodes(xin+1,:); % Binary integer form for m=1:size(grayCodes,1) map{m}=removeBlanks(num2str(grayCodes(m,:))); % String form end for m=1:numel(xin) grayCode_s(m,:)=map{xin(m)+1}; % String form end else fprintf(" Please specify the 2nd input argument signaling as one of 'ask'|'psk','qam'") end end function p=prob_error(SNRbdB,signaling,b,opt1,opt2) % Finds the symbol/bit error probability for given SNRbdB=Eb/(N0/2)[dB](Table 7.1) % Note that EbN0dB=SNRbdB-3. % Copyleft: Won Y. Yang, [email protected], CAU for academic use only if nargin<5, opt2='coherent'; end % opt2='coherent' or 'noncoherent' if nargin<4, opt1='BER'; end % opt1='SER' or 'BER' M=2^b; SNRb=10.^(SNRbdB/10); NSNR=length(SNRbdB); if signaling(1:3)=='ASK' % ASK (PAM) if lower(opt2(1))=='c' % ASK coherent p=2*(M-1)/M*Q(sqrt(3*b*SNRb/(M^2-1))); % Eq. (7.1.5) if lower(opt1(1))=='b', p = p/b; end else % ASK noncoherent if b==1, p=exp(-SNRb/4)/2; end %Eq.(7.1.22) end elseif signaling(1:3)=='FSK' tmp=M/2/(M-1); f5251_=@(x,SNRb,b)Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1); % Eq.(5.2.51) %f7208=inline('Q(-sqrt(2)*x-sqrt(b*SNRb)).^(2^b-1).*exp(-x.^2)','x','SNRb','b'); if lower(opt2(1))=='c' % FSK coherent if b==1 p=Q(sqrt(SNRb/2)); %Eq.(7.2.9) else Tol=1e-8; % tolerance used for 'quad' integration for i=1:NSNR p(i)=1-Gauss_Hermite(f5251_,10,SNRb(i),b)/sqrt(pi); %p(i)=1-adapt_smpsn('f7208',-20,20,Tol,SNRb(i),b)/sqrt(pi); % Eq.(7.2.8) %p(i)=1-quadl('f7208',-20,20,Tol,[],SNRb(i),b)/sqrt(pi); end end else % FSK noncoherent for i=1:NSNR p(i)=(M-1)/2*exp(-b*SNRb(i)/4); tmp1=M-1; for m=2:M-1 tmp1=-tmp1*(M-m)/m; p(i)=p(i)+tmp1/(m+1)*exp(-m*b*SNRb(i)/2/(m+1)); % Eq.(7.2.18) end end end if lower(opt1(1))=='b'&b>1, p=p*tmp; end % Eq.(7.2.19) elseif signaling(1:3)=='PSK' p=(1+(b>1))*Q(sqrt(b*SNRb)*sin(pi/M)); %Eq.(7.3.7) if lower(opt1(1))=='b'&b>1, p=p/b; end elseif signaling(1:3)=='DPS' %DPSK for i=1:NSNR EbN0=SNRb(i)/b; % From the MATLAB built-in function berawgn() tol=1e-6; % tol=max(min(1e-4/EbN0.^5,1e-4),eps); int_dpsk=@(y,EbN0,M,b)exp(-b*EbN0*(sin(pi/M))^2./(1+sqrt(1-(sin(pi/M)).^2).*cos(y))); %?????? p(i)=1/pi*quad(int_dpsk,1e-5,pi*(1-1/M),tol,[],EbN0,M,b); end if lower(opt1(1))=='b', p=p/b; end elseif signaling(1:3)=='QAM' L = 2^(ceil(b/2)); N = M/L; tmpL = 1-2*(L-1)/L*Q(sqrt(3*b/2/(L^2-1)*SNRb)); tmpN = 1-2*(N-1)/N*Q(sqrt(3*b/2/(N^2-1)*SNRb)); p = 1-tmpL.*tmpN; %Eq.(7.5.6) if lower(opt1(1))=='b'&b>1, p = p/b; end end end function newStr=removeBlanks(str) newStr = str(~isspace(str)); end

%cir06e01_1.m clear, clf Vm=450; Vs=Vm*exp(j*0); R1=40; R2=54; R3=40; L1=53.1e-3; L2=47.7e-3; C1=26.526e-6; C2=14.737e-6; f=60; w=2*pi*f; jw=j*w; % Phasor transform analysis Zjw=[R1+1/jw/C1 -1/jw/C1 -R1; -1/jw/C1 R2+R3+1/jw/C1 -R3; -R1 -R3 R1+R3+jw*L1]; Vjw=[Vs; 0; 0]; Ijw=Zjw\Vjw; IR3jw=Ijw(3)-Ijw(2) IR3m=abs(IR3jw), IR3p=angle(IR3jw) IR3p_deg=rad2deg(IR3p) tt=linspace(0,0.05,200); syms t iR3jw=IR3m*sin(w*t+IR3p) t=tt; iR3jwt=eval(iR3jw); plot(tt,iR3jwt), hold on % Laplace transform analysis syms s t Vs=Vm*w/(s^2+w^2); % Laplace transform of Vs Z=[R1+1/s/C1 -1/s/C1 -R1; -1/s/C1 R2+R3+1/s/C1 -R3; -R1 -R3 R1+R3+s*L1];; V=[Vs; 0; 0]; I=Z\V; IR3=I(3)-I(2); disp(IR3), pretty(IR3) iR3=ilaplace(IR3); disp(iR3'), pretty(iR3) iR3t=subs(iR3,t,tt); plot(tt,iR3t,'r') iR3_=ilaplace_my(IR3) iR3_=str2syms(iR3_); pretty(iR3_) iR3t_=subs(iR3_,t,tt); plot(tt,iR3t_,'g.')

%LPF_design_stubs_1_0.m % Determine the normalized LPF prototype element values: Rs=50; % Source (line) characteristic impedance N=3; type='C'; Rp=3; % for 3rd order Chebyshev with 3dB passband ripple gg=LPF_ladder_coef(N,type,Rp) % Determine the normalized LC values of the prototype LPF: Rs0=1; pT='T'; wc=1; % for a T-type prototype filter [LCRs,LCRs_,Gs]=LC_ladder(Rs0,gg,pT,wc); LCRs, LCRs_, pretty(simplify(Gs)) % Plot the input amplitude response |G(jω)|: fc=4e9; fN=linspace(0,pi/4,1000); % Normalized frequency range f=fN*2*pi*fc; % Real frequency range s=j*f/fc; Gw=eval(Gs); Gwmag=20*log10(max(abs(Gw),1e-5)); figure(1), clf plot(f,Gwmag-max(Gwmag)) % To make the maximum 0dB % Convert the L/C values into the impedances of short/open-circuited stubs: [Zs,SPsos]=LC2Z(LCRs,LCRs_) % Use Kuroda identity 1 to convert Z1 and Z5 into Ss’s. [Z0s(1),Zs(1),SPsos{1}]=Kuroda(Zs(1),SPsos{1}); [Z0s(3),Zs(3),SPsos{3}]=Kuroda(Zs(3),SPsos{3}); Z0s, Zs, SPsos % Denormalize the normalized characteristic impedances % by multiplying them with Rs=50Ω: Z0s=Z0s*Rs, Zs=Zs*Rs % Draw the designed microstrip filter. figure(2), clf draw_ustrip_fiter(Z0s,Zs), shg function [gg,a,b]=LPF_ladder_coef(N,type,Rp) % Input: % N = Filter order % type = 'B'/'L'/'C' for Butterworth-3dB/Linear-phase/Chebyshev-Rp[dB] % Output: % gg = Normalized LPF prototype LC element values % a, b = ??????????????????????????????????????????????????? % Ref[1]: "RF and Microwave Circuit and Component Design for Wireless Systems" % Kai Chang, Inder Bahl, and Vijay Nair (Wiley, 2002) % Ref[2]: "RF Circuit Design: Theory and Applications" % Reinhold Ludwig and Gene Bogdanov (Prentice Hall, 2009) if nargin<3, Rp=3; end % Rp is used ony for Chebyshev filters L=[2 1 0 0 0 0 0 0 0 0 0; 1.5774 0.4226 1 0 0 0 0 0 0 0 0; 1.255 0.5528 0.1922 1 0 0 0 0 0 0 0; 1.0598 0.5116 0.3181 0.1104 1 0 0 0 0 0 0; 0.9303 0.4577 0.3312 0.209 0.0718 1 0 0 0 0 0; 0.8377 0.4116 0.3158 0.2364 0.248 0.0505 1 0 0 0 0; 0.7677 0.3744 0.2944 0.2378 0.1778 0.1104 0.0375 1 0 0 0; 0.7125 0.3446 0.2735 0.2297 0.1867 0.1387 0.0855 0.0289 1 0 0; 0.6678 0.3203 0.2547 0.2184 0.1859 0.1506 0.1111 0.0682 0.023 1 0; 0.6305 0.3002 0.2384 0.2066 0.1808 0.1539 0.124 0.0911 0.0557 0.0187 1]; % Linear-Phase, Table 5.3 of Ref[2] % Data by courtesy of L. Weinberg and the Journal of the Franklin Institute. k=1:N; T=upper(type(1)); if T=='B' % Butterworth gg=[2*sin((2*k-1)*pi/2/N) 1]; % Eq.(6.8) of Ref[1] or Table 5-2 of Ref[2] % or <ecee.colorado.edu/~mathys/ecen2420/pdf/UsingFilterTables.pdf> elseif T=='L', gg=L(N,1:N+1); % Linear-phase elseif T=='C' % Table 5-4(a)/(b) for Chebyshev with Rp=3/0.5[dB], Ref[2] beta=log(coth(Rp/17.37)); gamma=sinh(beta/2/N); % Table 8.1 of Ref[1] a=sin((2*k-1)*pi/2/N); b=gamma^2+sin(k*pi/N).^2; % Eq.(6.11) of Ref[1] gg(1)=2*a(1)/gamma; % for k=2:N, gg(k)=4*a(k-1)*a(k)/b(k-1)/gg(k-1); end % Eq.(6.10) of Ref[1] if mod(N,2)==1, gg(N+1)=1; else gg(N+1)=coth(beta/4)^2; end end end function [LCRs,LCRs_,Gs]=LC_ladder(Rs,gg,pT,wc,band) % Input: % Rs = Source resistance % gg = Ladder coefficients generated using LPF_ladder_coef() % pT = 'p'/'T' or 1/2 for Pi/T-ladder % wc = wc/[wl wu] for band=('LPF'|'HPF')/('BPF'|'BSF') % band='LPF'/'BPF'/'BSF'/'HPF' % Output: % LCRs,LCRs_ = [C1 L2 C3 RL],{'C(v)','L(h)','C(v)' 'R(v)'} for Pi-ladder like Fig. 6.14.1(a) % LCRs,LCRs_ = [L1 C2 L3 RL],{'L(h)','C(v)','L(h)' 'R(v)'} for T-ladder like Fig. 6.14.1(b) % LCRs,LCRs_ = [L1 C1 L2 C2 L3 C3 RL],{'L-C(v)','L|C(h)','L-C(v)' 'R(v)'} for Pi-ladder like Fig. 6.14.4(a) % LCRs,LCRs_ = [L1 C1 L2 C2 L3 C3 RL],{'L|C(h)','L-C(v)','L|C(h)' 'R(v)'} for T-ladder like Fig. 6.14.4(b) % Gs = G(s): Transfer function % Copyleft: Won Y. Yang, [email protected], CAU for academic use only % LC Butterworth filter calculator: % <www.daycounter.com/Filters/LC-Butterworth-Filter/LC-Butterworth-Filter-Calculator.phtml> if nargin<5, band='LPF'; end if nargin<4, wc=1; end if length(wc)>1, w0=sqrt(wc(1)*wc(2)); wb=wc(2)-wc(1); end syms s N1=length(gg); N=N1-1; if isnumeric(pT), Kc=pT; else pT=lower(pT); Kc=(pT(1)=='p')+2*(pT(1)=='t'); % Set Kc=1/2 for pi/T-type end if mod(Kc+N,2)==1, RL=gg(N1)*Rs; % When gg(N1)=RL/Rs else RL=Rs/gg(N1); % When gg(N1)=Rs/RL end I=1; V=RL; % Set the current through RL to 1 to find the transfer ftn G(s) LCs=[]; LCs_=[]; % Set LCs and LCs_ as empty vectors. if Kc==1 % Pi-ladder for n=N:-1:1 if mod(n,2)==0, [LCsn,LCs_{n},V]=LC_horizontal(Rs,gg(n),wc,band,V,I); else [LCsn,LCs_{n},I]=LC_vertical(Rs,gg(n),wc,band,V,I); end LCs=[LCsn LCs]; end else % T-ladder for n=N:-1:1 if mod(n,2)==1, [LCsn,LCs_{n},V]=LC_horizontal(Rs,gg(n),wc,band,V,I); else [LCsn,LCs_{n},I]=LC_vertical(Rs,gg(n),wc,band,V,I); end LCs=[LCsn LCs]; end end LCRs=[LCs RL]; LCRs_=LCs_; LCRs_{N1}='R(v)'; Gs=RL/(V+Rs*I); % Transfer function %Gs=RL/V; % If Rs is not considered, i.e., Rs=0 end function [Z0,Z2,SPso]=Kuroda(Z1,SPso,Z0) % Series-to-Parallel of Z1(L) or Z1(short) for SPso='SL' or 'SS' % Series-to-Parallel of Z1(C) or Z1(open) for SPso='SC' or 'SO' % Parallel-to-Series of Z1(L) or Z1(short) for SPso='PL' or 'PS' % Parallel-to-Series of Z1(C) or Z1(open) for SPso='PC' or 'PO' % where Z0 is the characteristic impedance of a UE (unit element) % of electrical length bl=pi/4. % Input: % Z = Impedance (Reactance) of L or C % Z0 = Characteristic impedance of a UE if nargin<3|~isnumeric(Z0)|Z0==0, Z0=1; end if nargin<2, SPsp='SS'; end SPso=upper(SPso); if SPso(1)=='S' % Series-to-Parallel conversion if SPso(2)=='L'|SPso(2)=='S', SPso='Po'; n=1+Z0/Z1; Z2=-n*Z0; Z0=n*Z1; % Kuroda 2 else SPso='Po'; Z1=abs(Z1); n=1+Z1/Z0; Z0=n*Z0; Z2=[-n*Z1 n]; % Kuroda 4 end % The +/- sign of Z2 means L(short stub)/C(open stub), respectively. else % Parallel-to-Series conversion if SPso(2)=='L'|SPso(2)=='S', SPso='Ss'; n=1+Z0/Z1; Z0=Z0/n; Z2=[Z1/n 1/n]; % Kuroda 3 else SPso='Ss'; Z1=abs(Z1); n=1+Z1/Z0; Z2=Z0/n; Z0=Z1/n; % Kuroda 1 end end end function draw_ustrip_fiter(Z0s,Zs) Z0s1=nonzeros(Z0s); M=numel(Z0s1); Mh=floor(numel(Z0s)/2); Zs1=abs(Zs); N=numel(Zs); m0=0; n0=0; Ohm=char(937); for m=1:M plot([m-1 m],[1 1],'b', m-1,1,'bo', m,1,'bo'), hold on text(m-0.7,1.1,[num2str(Z0s1(m),'%.1f') Ohm]) plot([m-1 m],[0 0],'b', m-1,0,'bo', m,0,'bo') end for n=1:N plot([n-1 n-0.3],[1 0.7],'r', n-0.3,0.7,'ro') plot([n-1 n-0.3],[0 -0.3],'r', n-0.3,-0.3,'ro') text(n-0.9,0.5,[num2str(Zs1(n),'%.1f') Ohm]) end axis('equal') ax1=gca; % gca = get current axis ax1.YAxis.Visible='off'; % remove y-axis ax1.XAxis.Visible='off'; % remove x-axis end

Watch hkzbin.info/www/bejne/eXXUdYp4edqohtE for another design/implementation of RF microchip filter with stubs using MATLAB. The M-file listed below can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/173790-design-implementation-of-rf-microstrip-filter-with-stubs. %elec07e05.m clear, clf Rs=50; % Source (line) characteristic impedance N=3; type='C'; Rp=3; % for 3rd order Chebyshev with 3dB passband ripple gg=LPF_ladder_coef(N,type,Rp) Rs0=1; pT='T'; wc=1; % for a T-type prototype filter [LCRs,LCRs_,Gs]=LC_ladder(Rs0,gg,pT,wc); LCRs, LCRs_, pretty(simplify(Gs)) [Zs,SPsos]=LC2Z(LCRs,LCRs_) Zs_0=Zs; SPsos_0=SPsos; % To save the original values of Zs and SPsos [Z0s(1),Zs(1),SPsos{1}]=Kuroda(Zs(1),SPsos{1}); [Z0s(3),Zs(3),SPsos{3}]=Kuroda(Zs(3),SPsos{3}); Z0s, Zs, SPsos % Alternatively, using S2P_stubs() Zs=Zs_0; Z0s=zeros(size(Zs)); SPsos=SPsos_0; [Z0s_f,Zs_f,SPsos_f]=S2P_stubs(Z0s(1),Zs(1),SPsos(1)); [Z0s_b,Zs_b,SPsos_b]=S2P_stubs(Z0s(3),Zs(3),SPsos(3)); Z0s=[Z0s_f Z0s(2) Z0s_b(end:-1:1)], Zs=[Zs_f Zs(2) Zs_b(end:-1:1)] SPsos=[SPsos_f SPsos(2) SPsos_b(end:-1:1)] % Denormalize the normalized characteristic impedances % by multiplying them with Rs=50¥Ø. Z0s=Z0s*Rs, Zs=Zs*Rs % To plot the frequency response of the filter Z0strip_=[Zs(1) Z0s(1) Zs(2) Z0s(3) Zs(3)] SPsos_={SPsos{1} 'Tl' SPsos{2} 'Tl' SPsos{3}} [Z0strip,SPsos,LCRs,LCRs_]=ustrip_LPF(Rs,type,Rp,N,pT) l=1/8; fc=4e9; % Stub length and Cutoff frequency fN=linspace(0,pi/4,1000); % Normalized frequency range w=fN*2*pi*fc; % Real frequency range bl=2*pi*w*l/fc; % Range of bl (bl=2*pi*l at f=fc=4e9[Hz]) % To get the frequency response and input impedance ZL=Rs; [Gw,Zinw]=GZ_strip(Z0strip,SPsos,ZL,bl); Gw=Zinw./(Rs+Zinw).*Gw; % To get the overall frequency response considering Rs Gwmag=20*log10(max(abs(Gw),1e-5)); figure(1), clf plot(w,Gwmag-max(Gwmag)), hold on s=j*w/fc; Gw0=eval(Gs); Gwmag0=20*log10(max(abs(Gw0),1e-5)); plot(w,Gwmag0-max(Gwmag0),'r:') % To make the maximum 0dB legend('Microstrip','Lumped L/C') axis([w([1 end]) -50 5]) set(gcf,'color','white') figure(2), clf draw_ustrip_fiter(Z0s,Zs), shg function [gg,a,b]=LPF_ladder_coef(N,type,Rp) % Input: % N = Filter order % type = 'B'/'L'/'C' for Butterworth-3dB/Linear-phase/Chebyshev-Rp[dB] % Output: % gg = Normalized LPF prototype LC element values % a, b = ??????????????????????????????????????????????????? % Ref[1]: "RF and Microwave Circuit and Component Design for Wireless Systems" % Kai Chang, Inder Bahl, and Vijay Nair (Wiley, 2002) % Ref[2]: "RF Circuit Design: Theory and Applications" % Reinhold Ludwig and Gene Bogdanov (Prentice Hall, 2009) if nargin<3, Rp=3; end % Rp is used ony for Chebyshev filters L=[2 1 0 0 0 0 0 0 0 0 0; 1.5774 0.4226 1 0 0 0 0 0 0 0 0; 1.255 0.5528 0.1922 1 0 0 0 0 0 0 0; 1.0598 0.5116 0.3181 0.1104 1 0 0 0 0 0 0; 0.9303 0.4577 0.3312 0.209 0.0718 1 0 0 0 0 0; 0.8377 0.4116 0.3158 0.2364 0.248 0.0505 1 0 0 0 0; 0.7677 0.3744 0.2944 0.2378 0.1778 0.1104 0.0375 1 0 0 0; 0.7125 0.3446 0.2735 0.2297 0.1867 0.1387 0.0855 0.0289 1 0 0; 0.6678 0.3203 0.2547 0.2184 0.1859 0.1506 0.1111 0.0682 0.023 1 0; 0.6305 0.3002 0.2384 0.2066 0.1808 0.1539 0.124 0.0911 0.0557 0.0187 1]; % Linear-Phase, Table 5.3 of Ref[2] % Data by courtesy of L. Weinberg and the Journal of the Franklin Institute. k=1:N; T=upper(type(1)); if T=='B' % Butterworth gg=[2*sin((2*k-1)*pi/2/N) 1]; % Eq.(6.8) of Ref[1] or Table 5-2 of Ref[2] % or <ecee.colorado.edu/~mathys/ecen2420/pdf/UsingFilterTables.pdf> elseif T=='L', gg=L(N,1:N+1); % Linear-phase elseif T=='C' % Table 5-4(a)/(b) for Chebyshev with Rp=3/0.5[dB], Ref[2] beta=log(coth(Rp/17.37)); gamma=sinh(beta/2/N); % Table 8.1 of Ref[1] a=sin((2*k-1)*pi/2/N); b=gamma^2+sin(k*pi/N).^2; % Eq.(6.11) of Ref[1] gg(1)=2*a(1)/gamma; % for k=2:N, gg(k)=4*a(k-1)*a(k)/b(k-1)/gg(k-1); end % Eq.(6.10) of Ref[1] if mod(N,2)==1, gg(N+1)=1; else gg(N+1)=coth(beta/4)^2; end end end function [LCRs,LCRs_,Gs]=LC_ladder(Rs,gg,pT,wc,band) % Input: % Rs = Source resistance % gg = Ladder coefficients generated using LPF_ladder_coef() % pT = 'p'/'T' or 1/2 for Pi/T-ladder % wc = wc/[wl wu] for band=('LPF'|'HPF')/('BPF'|'BSF') % band='LPF'/'BPF'/'BSF'/'HPF' % Output: % LCRs,LCRs_ = [C1 L2 C3 RL],{'C(v)','L(h)','C(v)' 'R(v)'} for Pi-ladder like Fig. 6.14.1(a) % LCRs,LCRs_ = [L1 C2 L3 RL],{'L(h)','C(v)','L(h)' 'R(v)'} for T-ladder like Fig. 6.14.1(b) % LCRs,LCRs_ = [L1 C1 L2 C2 L3 C3 RL],{'L-C(v)','L|C(h)','L-C(v)' 'R(v)'} for Pi-ladder like Fig. 6.14.4(a) % LCRs,LCRs_ = [L1 C1 L2 C2 L3 C3 RL],{'L|C(h)','L-C(v)','L|C(h)' 'R(v)'} for T-ladder like Fig. 6.14.4(b) % Gs = G(s): Transfer function % Copyleft: Won Y. Yang, [email protected], CAU for academic use only % LC Butterworth filter calculator: % <www.daycounter.com/Filters/LC-Butterworth-Filter/LC-Butterworth-Filter-Calculator.phtml> if nargin<5, band='LPF'; end if nargin<4, wc=1; end if length(wc)>1, w0=sqrt(wc(1)*wc(2)); wb=wc(2)-wc(1); end syms s N1=length(gg); N=N1-1; if isnumeric(pT), Kc=pT; else pT=lower(pT); Kc=(pT(1)=='p')+2*(pT(1)=='t'); % Set Kc=1/2 for pi/T-type end if mod(Kc+N,2)==1, RL=gg(N1)*Rs; % When gg(N1)=RL/Rs else RL=Rs/gg(N1); % When gg(N1)=Rs/RL end I=1; V=RL; % Set the current through RL to 1 to find the transfer ftn G(s) LCs=[]; LCs_=[]; % Set LCs and LCs_ as empty vectors. if Kc==1 % Pi-ladder for n=N:-1:1 if mod(n,2)==0, [LCsn,LCs_{n},V]=LC_horizontal(Rs,gg(n),wc,band,V,I); else [LCsn,LCs_{n},I]=LC_vertical(Rs,gg(n),wc,band,V,I); end LCs=[LCsn LCs]; end else % T-ladder for n=N:-1:1 if mod(n,2)==1, [LCsn,LCs_{n},V]=LC_horizontal(Rs,gg(n),wc,band,V,I); else [LCsn,LCs_{n},I]=LC_vertical(Rs,gg(n),wc,band,V,I); end LCs=[LCsn LCs]; end end LCRs=[LCs RL]; LCRs_=LCs_; LCRs_{N1}='R(v)'; Gs=RL/(V+Rs*I); % Transfer function %Gs=RL/V; % If Rs is not considered, i.e., Rs=0 end function [Z0,Z2,SPso]=Kuroda(Z1,SPso,Z0) % Series-to-Parallel of Z1(L) or Z1(short) for SPso='SL' or 'SS' % Series-to-Parallel of Z1(C) or Z1(open) for SPso='SC' or 'SO' % Parallel-to-Series of Z1(L) or Z1(short) for SPso='PL' or 'PS' % Parallel-to-Series of Z1(C) or Z1(open) for SPso='PC' or 'PO' % where Z0 is the characteristic impedance of a UE (unit element) % of electrical length bl=pi/4. % Input: % Z = Impedance (Reactance) of L or C % Z0 = Characteristic impedance of a UE if nargin<3|~isnumeric(Z0)|Z0==0, Z0=1; end if nargin<2, SPsp='SS'; end SPso=upper(SPso); if SPso(1)=='S' % Series-to-Parallel conversion if SPso(2)=='L'|SPso(2)=='S', SPso='Po'; n=1+Z0/Z1; Z2=-n*Z0; Z0=n*Z1; % Kuroda 2 else SPso='Po'; Z1=abs(Z1); n=1+Z1/Z0; Z0=n*Z0; Z2=[-n*Z1 n]; % Kuroda 4 end % The +/- sign of Z2 means L(short stub)/C(open stub), respectively. else % Parallel-to-Series conversion if SPso(2)=='L'|SPso(2)=='S', SPso='Ss'; n=1+Z0/Z1; Z0=Z0/n; Z2=[Z1/n 1/n]; % Kuroda 3 else SPso='Ss'; Z1=abs(Z1); n=1+Z1/Z0; Z2=Z0/n; Z0=Z1/n; % Kuroda 1 end end end function [Gw,Zinw]=GZ_strip(Z0strip,SPsos,ZL,bl) % To get the frequency response Gw and input impedance Zinw % Input: % Z0strip = Impedances of stub, line, stub, line, ... % SPsos = e.g. {'Ss','Po',..} for Series-short, Parallel-open, .. % ZL = Load impedance % bl = Range of electrical length % Output: % Gw = Frequency response % Zinw = Input AC impedance Zl=@(ZL,Z0,bl)(ZL+j*Z0*tan(bl))./(Z0+j*ZL.*tan(bl))*Z0; % Impedance of a line of characteristic impedance Z0 and electrical length bl with load ZL Zls=@(Z0,bl)j*Z0*tan(bl); % Impedance of a short-circuited stub Zlo=@(Z0,bl)-j*Z0./tan(max([bl; eps*ones(1,numel(bl))])); % Impedance of an open-circuited stub Zlso=@(Z0,bl,so)(so=='s')*Zls(Z0,bl)+(so=='o')*Zlo(Z0,bl); Z_par=@(Z1,Z2)Z1.*Z2./(Z1+Z2); % Parallel(shunt) combination Z0strip=abs(Z0strip); Zinw=ZL; I=1; for m=numel(Z0strip):-1:1 if isempty(SPsos{m})|Z0strip(m)<eps, continue; end SP=upper(SPsos{m}(1)); so=lower(SPsos{m}(2)); if SP=='P' Zp=Z_par(Zlso(Z0strip(m),bl,so),Zinw); % Parallel combination I=I.*Zinw./Zp; Zinw=Zp; % Line current must be the sum of two shunt currents. elseif SP=='S', Zinw=Zinw+Zlso(Z0strip(m),bl,so); % Series combination elseif SP=='T', Zinw=Zl(Zinw,Z0strip(m),bl); % Serial line end end Gw=ZL./(Zinw.*I); % Ratio of the voltage across ZL to Input voltage end function draw_ustrip_fiter(Z0s,Zs) Z0s1=nonzeros(Z0s); M=numel(Z0s1); Mh=floor(numel(Z0s)/2); Zs1=abs(Zs); N=numel(Zs); m0=0; n0=0; Ohm=char(937); for m=1:M plot([m-1 m],[1 1],'b', m-1,1,'bo', m,1,'bo'), hold on text(m-0.7,1.1,[num2str(Z0s1(m),'%.1f') Ohm]) plot([m-1 m],[0 0],'b', m-1,0,'bo', m,0,'bo') end for n=1:N plot([n-1 n-0.3],[1 0.7],'r', n-0.3,0.7,'ro') plot([n-1 n-0.3],[0 -0.3],'r', n-0.3,-0.3,'ro') text(n-0.9,0.5,[num2str(Zs1(n),'%.1f') Ohm]) end axis('equal') ax1=gca; % gca = get current axis ax1.YAxis.Visible='off'; % remove y-axis ax1.XAxis.Visible='off'; % remove x-axis end

Watch kzbin.info/www/bejne/gn-3k6uuhKZ7g9U for another design/implementation of RF microchip filter with stubs using MATLAB. The m-file listed below can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/173795-design-implementation-of-rf-microstrip-filter. %LPF_design_stubs_2.m % Determine the normalized LPF prototype element values: wp=2*pi*1.6e9; ws=2*pi*3e9; Rp=0.5; As=15; [N,wc]=buttord(wp,ws,Rp,As,'s'); % Order and Cutoff freq of an analog Butterworth filter N, fc=wc/2/pi % Filter order and Cutoff (3dB) frequency[Hz] % Determine the normalized LPF prototype element values: Rs=50; % Source (line) characteristic impedance type='B'; % for a Butterworth LPF gg=LPF_ladder_coef(N,type,Rp) % Filter coefficients % Determine the normalized LC values of the prototype LPF: Rs0=1; pT='pi'; wc0=1; % for a pi-type prototype filter [LCRs,LCRs_,Gs]=LC_ladder(Rs0,gg,pT,wc0); LCRs, LCRs_, pretty(simplify(Gs)) % Plot the input amplitude response |G(jω)|: fN=linspace(0,pi/4,1000); % Normalized frequency range f=fN*2*pi*fc; % Real frequency range s=j*f/fc; Gw=eval(Gs); Gwmag=20*log10(max(abs(Gw),1e-5)); figure(1), clf plot(f,Gwmag-max(Gwmag)) % To make the maximum 0dB % Convert the L/C values into the impedances of short/open-circuited stubs: [Zs,SPsos]=LC2Z(LCRs,LCRs_) % Use Kuroda identity 1 to convert Z1 and Z5 into Ss’s. [Z0s(2),Zs(1),SPsos{1}]=Kuroda(Zs(1),SPsos{1}); [Z0s(4),Zs(5),SPsos{5}]=Kuroda(Zs(5),SPsos{5}); Z0s, Zs, SPsos % Insert UEs of Z0=1 between them and TL as shown in Fig.(d). Z0s(1)=1; Z0s(5)=1; % Use Kuroda identity 2 to convert the Ss’s Z1, Z2, Z4, & Z5 into Po’s. for m=[1 2 4 5] [Z0s(m),Zs(m),SPsos{m}]=Kuroda(Zs(m),SPsos{m},Z0s(m)); end Z0s, Zs, SPsos % Denormalize the normalized characteristic impedances % by multiplying them with Rs=50Ω: Z0s=Z0s*Rs, Zs=Zs*Rs % Draw the designed microstrip filter. figure(2), clf draw_ustrip_fiter(Z0s,Zs), shg function [gg,a,b]=LPF_ladder_coef(N,type,Rp) % Input: % N = Filter order % type = 'B'/'L'/'C' for Butterworth-3dB/Linear-phase/Chebyshev-Rp[dB] % Output: % gg = Normalized LPF prototype LC element values % a, b = ??????????????????????????????????????????????????? % Ref[1]: "RF and Microwave Circuit and Component Design for Wireless Systems" % Kai Chang, Inder Bahl, and Vijay Nair (Wiley, 2002) % Ref[2]: "RF Circuit Design: Theory and Applications" % Reinhold Ludwig and Gene Bogdanov (Prentice Hall, 2009) if nargin<3, Rp=3; end % Rp is used ony for Chebyshev filters L=[2 1 0 0 0 0 0 0 0 0 0; 1.5774 0.4226 1 0 0 0 0 0 0 0 0; 1.255 0.5528 0.1922 1 0 0 0 0 0 0 0; 1.0598 0.5116 0.3181 0.1104 1 0 0 0 0 0 0; 0.9303 0.4577 0.3312 0.209 0.0718 1 0 0 0 0 0; 0.8377 0.4116 0.3158 0.2364 0.248 0.0505 1 0 0 0 0; 0.7677 0.3744 0.2944 0.2378 0.1778 0.1104 0.0375 1 0 0 0; 0.7125 0.3446 0.2735 0.2297 0.1867 0.1387 0.0855 0.0289 1 0 0; 0.6678 0.3203 0.2547 0.2184 0.1859 0.1506 0.1111 0.0682 0.023 1 0; 0.6305 0.3002 0.2384 0.2066 0.1808 0.1539 0.124 0.0911 0.0557 0.0187 1]; % Linear-Phase, Table 5.3 of Ref[2] % Data by courtesy of L. Weinberg and the Journal of the Franklin Institute. k=1:N; T=upper(type(1)); if T=='B' % Butterworth gg=[2*sin((2*k-1)*pi/2/N) 1]; % Eq.(6.8) of Ref[1] or Table 5-2 of Ref[2] % or <ecee.colorado.edu/~mathys/ecen2420/pdf/UsingFilterTables.pdf> elseif T=='L', gg=L(N,1:N+1); % Linear-phase elseif T=='C' % Table 5-4(a)/(b) for Chebyshev with Rp=3/0.5[dB], Ref[2] beta=log(coth(Rp/17.37)); gamma=sinh(beta/2/N); % Table 8.1 of Ref[1] a=sin((2*k-1)*pi/2/N); b=gamma^2+sin(k*pi/N).^2; % Eq.(6.11) of Ref[1] gg(1)=2*a(1)/gamma; % for k=2:N, gg(k)=4*a(k-1)*a(k)/b(k-1)/gg(k-1); end % Eq.(6.10) of Ref[1] if mod(N,2)==1, gg(N+1)=1; else gg(N+1)=coth(beta/4)^2; end end end function [LCRs,LCRs_,Gs]=LC_ladder(Rs,gg,pT,wc,band) % Input: % Rs = Source resistance % gg = Ladder coefficients generated using LPF_ladder_coef() % pT = 'p'/'T' or 1/2 for Pi/T-ladder % wc = wc/[wl wu] for band=('LPF'|'HPF')/('BPF'|'BSF') % band='LPF'/'BPF'/'BSF'/'HPF' % Output: % LCRs,LCRs_ = [C1 L2 C3 RL],{'C(v)','L(h)','C(v)' 'R(v)'} for Pi-ladder like Fig. 6.14.1(a) % LCRs,LCRs_ = [L1 C2 L3 RL],{'L(h)','C(v)','L(h)' 'R(v)'} for T-ladder like Fig. 6.14.1(b) % LCRs,LCRs_ = [L1 C1 L2 C2 L3 C3 RL],{'L-C(v)','L|C(h)','L-C(v)' 'R(v)'} for Pi-ladder like Fig. 6.14.4(a) % LCRs,LCRs_ = [L1 C1 L2 C2 L3 C3 RL],{'L|C(h)','L-C(v)','L|C(h)' 'R(v)'} for T-ladder like Fig. 6.14.4(b) % Gs = G(s): Transfer function % Copyleft: Won Y. Yang, [email protected], CAU for academic use only % LC Butterworth filter calculator: % <www.daycounter.com/Filters/LC-Butterworth-Filter/LC-Butterworth-Filter-Calculator.phtml> if nargin<5, band='LPF'; end if nargin<4, wc=1; end if length(wc)>1, w0=sqrt(wc(1)*wc(2)); wb=wc(2)-wc(1); end syms s N1=length(gg); N=N1-1; if isnumeric(pT), Kc=pT; else pT=lower(pT); Kc=(pT(1)=='p')+2*(pT(1)=='t'); % Set Kc=1/2 for pi/T-type end if mod(Kc+N,2)==1, RL=gg(N1)*Rs; % When gg(N1)=RL/Rs else RL=Rs/gg(N1); % When gg(N1)=Rs/RL end I=1; V=RL; % Set the current through RL to 1 to find the transfer ftn G(s) LCs=[]; LCs_=[]; % Set LCs and LCs_ as empty vectors. if Kc==1 % Pi-ladder for n=N:-1:1 if mod(n,2)==0, [LCsn,LCs_{n},V]=LC_horizontal(Rs,gg(n),wc,band,V,I); else [LCsn,LCs_{n},I]=LC_vertical(Rs,gg(n),wc,band,V,I); end LCs=[LCsn LCs]; end else % T-ladder for n=N:-1:1 if mod(n,2)==1, [LCsn,LCs_{n},V]=LC_horizontal(Rs,gg(n),wc,band,V,I); else [LCsn,LCs_{n},I]=LC_vertical(Rs,gg(n),wc,band,V,I); end LCs=[LCsn LCs]; end end LCRs=[LCs RL]; LCRs_=LCs_; LCRs_{N1}='R(v)'; Gs=RL/(V+Rs*I); % Transfer function %Gs=RL/V; % If Rs is not considered, i.e., Rs=0 end function [Zs,SPsos]=LC2Z(LCs,LCs_) % To convert the L/C's into short/open-circuited stubs for m=1:numel(LCs) if LCs_{m}(1)=='L', Zs(m)=LCs(m); if LCs_{m}(3)=='h',SPsos{m}='Ss'; else SPsos{m}='Ps'; end elseif LCs_{m}(1)=='C', Zs(m)=-1/LCs(m); if LCs_{m}(3)=='h',SPsos{m}='So'; else SPsos{m}='Po'; end else break; end end end function [Z0,Z2,SPso]=Kuroda(Z1,SPso,Z0) % Series-to-Parallel of Z1(L) or Z1(short) for SPso='SL' or 'SS' % Series-to-Parallel of Z1(C) or Z1(open) for SPso='SC' or 'SO' % Parallel-to-Series of Z1(L) or Z1(short) for SPso='PL' or 'PS' % Parallel-to-Series of Z1(C) or Z1(open) for SPso='PC' or 'PO' % where Z0 is the characteristic impedance of a UE (unit element) % of electrical length bl=pi/4. % Input: % Z = Impedance (Reactance) of L or C % Z0 = Characteristic impedance of a UE if nargin<3|~isnumeric(Z0)|Z0==0, Z0=1; end if nargin<2, SPsp='SS'; end SPso=upper(SPso); if SPso(1)=='S' % Series-to-Parallel conversion if SPso(2)=='L'|SPso(2)=='S', SPso='Po'; n=1+Z0/Z1; Z2=-n*Z0; Z0=n*Z1; % Kuroda 2 else SPso='Po'; Z1=abs(Z1); n=1+Z1/Z0; Z0=n*Z0; Z2=[-n*Z1 n]; % Kuroda 4 end % The +/- sign of Z2 means L(short stub)/C(open stub), respectively. else % Parallel-to-Series conversion if SPso(2)=='L'|SPso(2)=='S', SPso='Ss'; n=1+Z0/Z1; Z0=Z0/n; Z2=[Z1/n 1/n]; % Kuroda 3 else SPso='Ss'; Z1=abs(Z1); n=1+Z1/Z0; Z2=Z0/n; Z0=Z1/n; % Kuroda 1 end end end function draw_ustrip_fiter(Z0s,Zs) Z0s1=nonzeros(Z0s); M=numel(Z0s1); Mh=floor(numel(Z0s)/2); Zs1=abs(Zs); N=numel(Zs); m0=0; n0=0; Ohm=char(937); for m=1:M plot([m-1 m],[1 1],'b', m-1,1,'bo', m,1,'bo'), hold on text(m-0.7,1.1,[num2str(Z0s1(m),'%.1f') Ohm]) plot([m-1 m],[0 0],'b', m-1,0,'bo', m,0,'bo') end for n=1:N plot([n-1 n-0.3],[1 0.7],'r', n-0.3,0.7,'ro') plot([n-1 n-0.3],[0 -0.3],'r', n-0.3,-0.3,'ro') text(n-0.9,0.5,[num2str(Zs1(n),'%.1f') Ohm]) end axis('equal') ax1=gca; % gca = get current axis ax1.YAxis.Visible='off'; % remove y-axis ax1.XAxis.Visible='off'; % remove x-axis end

%em07e01.m f=@(x,y)0; g=@(x,y)0; % Eq.(E7.1.1) w.r.t. Eq.(7.1.1) % (Rectangular) Domain and BC types x0=0; xf=pi; y0=0; yf=pi; D=[x0 xf y0 yf]; % BCs bx0=@(y)0; bxf=@(y)0; by0=@(x)sin(2*x); byf=@(x)0; % Eq.(E7.1.2a,b,c,d) % Numerical approach % Numbers of grid points along x-/y-axes, Error tolerance, Mx=50; My=50; tol=1e-5; MaxIter=600; [un,xx,yy]=pde_poisson_D(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter); % Analytical approach K=10; [ua,A]=pde_Laplace_ana(xf,yf,bx0,bxf,by0,byf,Mx,My,K); A, Discrepancy=norm(un-ua)/norm(ua) mesh(xx,yy,un) % Plot the solution graph xlabel('x'), ylabel('y') function [u,A,xx,yy]=pde_Laplace_ana(xf,yf,bx0,bxf,by0,byf,Mx,My,N) % To solve u_xx +u_yy = 0 % over the region D=[0,xf,0,yf]={(x,y)|0<=x<=xf, 0<=y<=yf} % with the boundary conditions: % u(0,y) = bx0(y), u(xf,y) = bxf(y) % u(x,0) = by0(x), u(x,yf) = byf(x) % Mx = Number of subintervals along x-axis % My = Number of subintervals along y-axis % N: Maximum number of terms if nargin<9, N=10; end if nargin<8, My=50; end if nargin<7, Mx=50; end x0=0; dx=(xf-x0)/Mx; xx=x0+[0:Mx]*dx; y0=0; dy=(yf-y0)/My; yy=y0+[0:My]'*dy; [X,Y]=meshgrid(xx,yy); u=zeros(size(X)); syms x y for n=1:N npixf = n*pi/xf; npiyf = n*pi/yf; lam(n)=yf*sinh(npiyf*xf); mu(n)=xf*sinh(npixf*yf); A(1,n)=2/lam(n)*eval(int(bx0*sin(npiyf*y),0,yf)); A(2,n)=2/lam(n)*eval(int(bxf*sin(npiyf*y),0,yf)); A(3,n)=2/mu(n)*eval(int(by0*sin(npixf*x),0,xf)); A(4,n)=2/mu(n)*eval(int(byf*sin(npixf*x),0,xf)); u = u + A(1,n)*sin(npiyf*Y).*sinh(npiyf*(xf-X)) + A(2,n)*sin(npiyf*Y).*sinh(npiyf*X) ... + A(3,n)*sin(npixf*X).*sinh(npixf*(yf-Y)) + A(4,n)*sin(npixf*X).*sinh(npixf*Y); end end function [u,x,y,itr]=pde_poisson_D(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter) % To solve u_xx +u_yy +g(x,y)u = f(x,y) % over the region D=[x0,xf,y0,yf]={(x,y)|x0<=x<=xf, y0<=y<=yf} % subject to the Dirichlet BCs (boundary conditions): % u(x0,y) = bx0(y), u(xf,y) = bxf(y) % u(x,y0) = by0(x), u(x,yf) = byf(x) % D = [x0 xf y0 yf] : Domain with Dirichlet BCs % Mx = Number of subintervals along x-axis % My = Number of subintervals along y-axis % tol : Error tolerance % MaxIter: Maximum number of iterations x0=D(1); xf=D(2); y0=D(3); yf=D(4); dx=(xf-x0)/Mx; x=x0+[0:Mx]*dx; dy=(yf-y0)/My; y=y0+[0:My]'*dy; Mx1=Mx+1; My1=My+1; % Dirichlet BCs (7.1.5b) for m=1:My1, u(m,[1 Mx1])=[bx0(y(m)) bxf(y(m))]; end % Left/Right side for n=1:Mx1, u([1 My1],n)=[by0(x(n)); byf(x(n))]; end % Bottom/Top % Initialize as the average of boundary values sum_of_bv=sum(sum([u(2:My,[1 Mx1]) u([1 My1],2:Mx)'])); u(2:My,2:Mx)=sum_of_bv/(2*(Mx+My-2)); for i=1:My for j=1:Mx, F(i,j)=f(x(j),y(i)); G(i,j)=g(x(j),y(i)); end end dx2=dx^2; dy2=dy^2; dxy2=2*(dx2+dy2); rx=dx2/dxy2; ry=dy2/dxy2; rxy=rx*dy2; for itr=1:MaxIter for j=2:Mx for i=2:My u(i,j)=ry*(u(i,j+1)+u(i,j-1)) +rx*(u(i+1,j)+u(i-1,j))... +rxy*(G(i,j)*u(i,j)-F(i,j)); % Eq.(7.1.5a) end end if itr>1 & max(max(abs(u-u0)))<tol, break; end u0=u; end end % Source: Engineering Mathematics with MATLAB by Won Y. Yang et al. (books.google.com) www.google.co.kr/books/edition/Engineering_Mathematics_with_MATLAB/R2OuEAAAQBAJ?hl=ko&gbpv=1&dq=Source:+Engineering+Mathematics+with+MATLAB++by+Won+Y.+Yang+et+al.+(books.google.com)&printsec=frontcover

%em07e03.m a=1e-2; it0=@(x,y)0; % IC (E7.3.2a) bxyt=@(x,y,t)exp(y)*cos(x)-exp(x)*cos(y); % BC (E7.3.2b) D=[0 4 0 4]; % Solution region t0=0; tf=5000; Mx=40; My=40; N=50; [u,x,y]=pde_heat2_ADI(a,D,tf,it0,bxyt,Mx,My,N); mesh(x,y,u) xlabel('x'), ylabel('y') function [u,x,y,t]=pde_heat2_ADI(a,D,tf,it0,bxyt,Mx,My,N) %To solve a^2*(u_xx+u_yy)=u_t for D(1)<=x<=D(2), D(3)<=y<=D(4), 0<=t<=tf % IC: u(x,y,0) = it0(x,y) % BC: u(x,y,t) = bxyt(x,y,t) for (x,y)cB % Mx/My = Number of subintervals along x/y-axis % N = Number of subintervals along t-axis dx=(D(2)-D(1))/Mx; x=D(1)+[0:Mx]*dx; dy=(D(4)-D(3))/My; y=D(3)+[0:My]'*dy; dt=tf/N; t=[0:N]*dt; % Initialization for j=1:Mx+1, for i=1:My+1, u(i,j)=it0(x(j),y(i)); end, end rx=a^2*dt/(dx*dx); rx1=1+2*rx; rx2=1-2*rx; ry=a^2*dt/(dy*dy); ry1=1+2*ry; ry2=1-2*ry; for j=1:Mx-1 % Eq.(7.2.38a) Ay(j,j)=ry1; if j>1, Ay(j-1,j)= -ry; Ay(j,j-1)= -ry; end end for i=1:My-1 % Eq.(7.2.38b) Ax(i,i)=rx1; if i>1, Ax(i-1,i)= -rx; Ax(i,i-1)= -rx; end end for k=1:N u_1=u; t=k*dt; for i=1:My+1, u(i,[1 Mx+1])=feval(bxyt,x([1 Mx+1]),y(i),t); end % BCs for j=1:Mx+1, u([1 My+1],j)=feval(bxyt,x(j),y([1;My+1]),t); end % BCs if mod(k,2)==0 % for k even for i=2:My jj=2:Mx; bx=[ry*u(i,1) zeros(1,Mx-3) ry*u(i,My+1)] ... +rx*(u_1(i-1,jj)+u_1(i+1,jj)) +rx2*u_1(i,jj); u(i,jj)= trid(Ay,bx')'; % Eq.(7.2.38a) end else % for k odd for j=2:Mx ii=2:My; by=[rx*u(1,j); zeros(My-3,1); rx*u(Mx+1,j)] ... +ry*(u_1(ii,j-1)+u_1(ii,j+1)) +ry2*u_1(ii,j); u(ii,j)= trid(Ax,by); % Eq.(7.2.38b) end end end

%inv_pendulum.m Mc=5; Mp=1; L=2; d=1; g=9.8; xr=[1; 0; 0; 0]; % Reference state A=[0 0 1 0; 0 0 0 1; 0 Mp*g/Mc -d/Mc 0; 0 (Mc+Mp)*g/Mc/L -d/Mc/L 0] B=[0; 0; 1/Mc; 1/Mc/L]; C=[1 0 0 0; 0 1 0 0]; Q=C'*C; %[1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0]; Q=1000*Q; % Increase these values to make the system response faster % possibly at the cost of larger force R=10; % Increase this value to reduce the input force K=lqr(A,B,Q,R); u=@(x)-K*(x-xr); % Control input (force) function x0=[0; 0.7; 0; 0]; % Initial vector t0=0; tf=10; tspan=[t0 tf]; % Time span of 10 sec % To solve the ODE [tt,xx]=ode45(@(t,x)dxdt(t,x,Mc,Mp,L,d,g,u(x)),tspan,x0); Nt=numel(tt); animation=1; if animation<1 ft=-K*(xx.'-repmat(xr,1,Nt)); %ft_=u(xx); subplot(3,1,1) plot(tt,xx(:,1),'b', tt([1 end]),xr(1)*[1 1],'m:') ylabel('x[m]') subplot(3,1,2) plot(tt,rad2deg(xx(:,2)),'r', tt([1 end]),rad2deg(xr(2))*[1 1],'m:') ylabel('theta[deg]') subplot(3,1,3), plot(tt,ft,'g') ylabel('force f(t)') else % Animation for n=1:Nt t=tt(n); tp=t; xn=xx(n,1); thn=xx(n,2); if n>1, pause(t-tp); cla; end draw_inv_pendulum(xn,thn,L) axis('equal'), axis([-12 8 -1 5]) end end function dx=dxdt(t,x,Mc,Mp,L,d,g,f) % State equation for an inverted pendulum %global Mc Mp L d g K xr th=x(2); w=x(4); sth=sin(th); cth=cos(th); A0=[Mc+Mp -Mp*L*cth; -cth L]; b0=[f-Mp*L*sth*w^2-d*x(3); g*sth]; %f=-K*(x-xr); dx34 = A0\b0; % Eq.(2) dx = [x(3); x(4); dx34(1); dx34(2)]; end function draw_inv_pendulum(x,th,L,w,h,color) %x : the x-coordinate of the cart %th: the angle[rad] of the pendulum cart from upright %h : the height of the cart %L : the length of the pendulum %w : the width of the cart %h : the height of the cart if nargin<6, color='k'; end if nargin<5, h=1.2*L; end if nargin<4, w=1.5*h; end xs1=x+[-w/2 -w/2 w/2 w/2 -w/2]; ys1=h/10*[1 10 10 1 1]; plot(xs1,ys1,color); hold on circle(x-w/3,h/10,h/10); circle(x+w/3,h/10,h/10); circle(x,h,h/20,-pi,0); px=x-L*sin(th); py=h+L*cos(th); plot([x px],[h py],'r') plot(px,py,'ro','MarkerSize',10) axis('equal') end

Thank you.

Here is the MATLAB code for solving the heat equation as shown above. %em07e02.m clear, clf a=1; % The parameter of Eq.(E7.2.1) it0=@(x)sin(pi*x); % Initial condition bx0=@(t)0; bxf=@(t)0; % Boundary condition xf=1; M=25; tf=0.1; N=120; % r=0.625 [u_exp,x,t]=pde_heat_exp(a,xf,tf,it0,bx0,bxf,M,N); [u_imp,x,t]=pde_heat_imp(a,xf,tf,it0,bx0,bxf,M,N); % converge unconditionally [u_CN,x,t]=pde_heat_CN(a,xf,tf,it0,bx0,bxf,M,N); % converge unconditionally subplot(131), mesh(t,x,u_exp), set(gca,'fontsize',9) title('Numerical solution using pde heat exp()') xlabel('t[sec]'), ylabel('x[m]') subplot(132), mesh(t,x,u_CN), set(gca,'fontsize',9) title('Numerical solution using pde heat CN()') xlabel('t[sec]'), ylabel('x[m]') m=0; sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); %u_pdepe = sol(:,:,1); % Extract the first solution component as u. u_pdepe = sol; % Extract the first solution component as u. % Analytical solution uo=@(x,t)sin(pi*x)*exp(-pi*pi*t); Uo= uo(x,t); aUo=abs(Uo)+eps; % Values of true analytical solution Nterm=10; [u_an,A,tt,xx]=pde_heat_ana(a,xf,tf,it0,M,N,Nterm); A %size(u_an), size(Uo) MN=M*N; err_an=norm((u_an-Uo)./aUo)/MN err_exp= norm((u_exp-Uo)./aUo)/MN err_imp= norm((u_imp-Uo)./aUo)/MN err_CN=norm((u_CN-Uo)./aUo)/MN err_pde=norm((u_pdepe.'-Uo)./aUo)/MN subplot(133), mesh(t,x,u_pdepe.') % surf(t,x,u) title('Numerical solution using pdepe()') set(gca,'fontsize',9), xlabel('t[sec]'), ylabel('x[m]') set(gcf,'color','white'), shg % -------------------------------------------------------------- function [c,f,s]=pdex1pde(x,t,u,DuDx) c=1; %1/a; f=DuDx; s=0; end function u0 = pdex1ic(x) u0=sin(pi*x); % Initial condition end function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl=ul; ql=0; % Left boundary condition pr=ur; qr=0; % Right boundary condition end function [u,x,t]=pde_heat_CN(a,xf,tf,it0,bx0,bxf,M,N) % To solve a^2*u_xx = u_t for 0<=x<=xf, 0<=t<=tf % Initial Condition: u(x,0) = it0(x) % Boundary Condition: u(0,t) = bx0(t), u(xf,t)=bxf(t) % M = Number of subintervals along x-axis % N = Number of subintervals along t-axis dx= xf/M; x= [0:M]'*dx; dt= tf/N; t= [0:N]*dt; for i=1:M+1, u(i,1)= it0(x(i)); end for n=1:N+1, u([1 M+1],n)= [bx0(t(n)); bxf(t(n))]; end %if nargin<9, Lmda=1; end r= a^2*dt/dx/dx r1=2*(1+r); r2=2*(1-r); for i=1:M-1 A(i,i)= r1; if i>1, A(i-1,i)= -r; A(i,i-1)= -r; end end for k=2:N+1 b= [r*u(1,k); zeros(M-3,1); r*u(M+1,k)] ... +r*(u(1:M-1,k-1)+u(3:M+1,k-1)) +r2*u(2:M,k-1); u(2:M,k)= trid(A,b); end end function [u,A,tt,xx]=pde_heat_ana(a,xf,tf,it0,M,N,Nterm) % To solve a^2*u_xx = u_t % over the region D=[0,tf,0,xf]={(t,x)|0<=t<=tf, 0<=x<=xf} % with the boundary/initial conditions: % u(t,0) = bx0(t), u(t,xf) = bxf(t), u(0,x) = it0(x) % M = Number of subintervals along x-axis % N = Number of subintervals along t-axis % Nterm: Maximum number of terms if nargin<8, Nterm=10; end if nargin<7, N=50; end if nargin<6, M=50; end x0=0; dx=(xf-x0)/M; xx= x0+[0:M]*dx; t0=0; dt=(tf-t0)/N; tt= t0+[0:N]'*dt; [T,X]=meshgrid(tt,xx); u=zeros(size(X)); syms x % x = sym('x') for n=1:Nterm npixf = n*pi/xf; A(n)=2/xf*eval(int(it0*sin(npixf*x),0,xf)); % Eq.(7.2.16) u = u + A(n)*sin(npixf*X).*exp(-a^2*npixf^2*T); % Eq.(7.2.14) end end function [u,x,t]=pde_heat_exp(a,xf,tf,it0,bx0,bxf,M,N) % To solve a^2*u_xx = u_t for 0<=x<=xf, 0<=t<=tf % Initial Condition: u(x,0) = it0(x) % Boundary Condition: u(0,t) = bx0(t), u(xf,t)=bxf(t) % M = Number of subintervals along x-axis % N = Number of subintervals along t-axis dx= xf/M; x= [0:M]'*dx; dt= tf/N; t= [0:N]*dt; for i=1:M+1, u(i,1)= it0(x(i)); end for n=1:N+1, u([1 M+1],n)= [bx0(t(n)); bxf(t(n))]; end r= a^2*dt/dx^2, r1= 1-2*r; for k=1:N for i=2:M u(i,k+1)= r*(u(i-1,k)+u(i+1,k)) +r1*u(i,k); end end end function [u,x,t]=pde_heat_imp(a,xf,tf,it0,bx0,bxf,M,N) % To solve a^2*u_xx = u_t for 0<=x<=xf, 0<=t<=tf % Initial Condition: u(x,0) = it0(x) % Boundary Condition: u(0,t) = bx0(t), u(xf,t)=bxf(t) % M = Number of subintervals along x-axis % N = Number of subintervals along t-axis dx= xf/M; x= [0:M]'*dx; dt= tf/N; t= [0:N]*dt; for i=1:M+1, u(i,1)= it0(x(i)); end for n=1:N+1, u([1 M+1],n)= [bx0(t(n)); bxf(t(n))]; end r= a^2*dt/dx^2; r2= 1+2*r; for i=1:M-1 A(i,i)= r2; if i>1, A(i-1,i)= -r; A(i,i-1)= -r; end %Eq.(9.2-9) end for k=2:N+1 b= [r*u(1,k); zeros(M-3,1); r*u(M+1,k)] +u(2:M,k-1); %Eq.(9.2-9) u(2:M,k)= trid(A,b); end end function x=trid(A,b) % Solve tridiagonal system of equations by Upper Triangularization % Applied Numerical Methods Using MATLAB (Wiley) by Won Y. Yang et. al N= size(A,2); row=0; if size(b,1)==1, b=b(:); row=1; end for m =2:N tmp =A(m,m-1)/A(m-1,m-1); A(m,m) =A(m,m) -A(m-1,m)*tmp; A(m,m-1) =0; b(m,:) =b(m,:) -b(m-1,:)*tmp; end % Back Substitution x(N,:) =b(N,:)/A(N,N); for m =N-1: -1: 1 x(m,:) =(b(m,:) -A(m,m+1)*x(m+1))/A(m,m); end if row==1, x=x'; end end end

%do_MBK_h.m clear figure(1), clf; figure(2), clf global M K B L g M=5; K=10; B=2; L=5; %g=9.8; %M=5; K=2; B=1; L=5; %g=9.8; t0=0; tf=20; x0=[8; 0]; [tt,xx]=ode45(@dxdt,[t0 tf],x0); animation=1; if animation<1 plot(tt,f(tt),'r', tt,xx(:,1),'b') else % Animatiom for n=1:numel(tt) t=tt(n); tp=t; figure(1), cla draw_MBK_h(4,1,xx(n,1),'b') axis('equal'), axis([-1 10 -3 3]) figure(2) plot(t,f(t),'r.', t,xx(n,1),'b.') axis([0 tf -3 10]), hold on if n>1, pause(t-tp); end end end shg function dx=dxdt(t,x) % State equation for an MKB system global M K B L dx=[x(2); (-K*(x(1)-L)-B*x(2)+f(t))/M]; end function f=f(t) % Force function f=10*sin(t); end function draw_MBK_h(n,w,x,color) %n: the # of spring windings %w: the width of each object %x: displacement of Mass if nargin<5, color='k'; end p1=[0 -w]; p2=[x -w]; p2_=[x-w -w]; xm=(p1(1)+p2(1))/2; ym=0 %Mass xM= p2(1)+w*[-1 1 1 -1 -1]; yM= ym+w*1.5*[-1 -1 1 1 -1]; plot([0 0],[min(yM) max(yM)],'k','LineWidth',2), hold on plot(xM,yM,color) shg %Spring spring_h(n,p1,p2_,w,color) %Damper damper_h(ym+w,p1(1),p2_(1),w,color) end function spring_h(n,p1,p2,w,color) %draw a horizontal spring of n windings, width w from p1 to p2 if nargin<5, color='k'; end L=w; % Length of spring end c= (p2(2)-p1(2))/2; d= (p2(1)-p1(1)-L)/2; f= (p2(2)+p1(2))/2; g= (p2(1)+p1(1))/2; x= -1:0.01:1; t= (x+1)*pi*(n+0.5); y=0.5*w*sin(t); x= x+0.15*(1-cos(t)); a=x(1); b=x(length(y)); x= 2*(x-a)/(b-a)-1; xxS= d*x - c*y +g; yyS= y + f; yyS1= [f f]; xxS1= xxS(end)+[0 L/2]; xxS2= xxS(1)-[0 L/2]; plot(xxS,yyS,color, xxS1,yyS1,color, xxS2,yyS1,color) end function damper_h(ym,x1,x2,w,color) %draws a damper in (ym-0.5 ym+0.5 x1 x2) if nargin<5, color='k'; end xm=(x1+x2)/2; dx=x2-x1; xD1= [x2 xm xm xm]; yD1= ym+w*0.2*[0 0 -1 1]; xD2= xm+dx/8*[1 -1 -1 1]; yD2= ym+w*0.3*[-1 -1 1 1]; xD3= [x1 xm-dx/8]; yD3= ym+[0 0]; plot(xD1,yD1,color, xD2,yD2,color, xD3,yD3,color) end

Do you have github repository of this code?

Thank you a lot. That is very useful for me to solve the electromagnetic problem

%do_MKB2.m clear, clf global M1 M2 K1 K2 B1 B2 M1=400; M2=40; K1=4; K2=1; B1=1000; B2=500; %M1=450; M2=40; K1=17500; K2=150000; B1=1500; B2=0; %M1=290; M2=60; K1=16200; K2=191000; B1=1000; B2=0; %f=@(t)0; % Active force %u=@(t)0.3*sin(4*t); du=@(t)1.2*cos(4*t); % Road disturbance and its derivative t0=0; tf=5; tspan=[t0 tf]; x0=[22; 10; 0; 0]; x0=[1.5; 0.5; 0; 0]; [tt,xx]=ode45(@dxdt,tspan,x0); Ms=[M1 M2]; Ks=[K1 K2]; Bs=[B1 B2]; animation = 1; w=0.2; if animation<1 plot(tt,u(tt), tt,xx(:,1:2)) legend('u(t)','y_1 (t)','y_2 (t)') else for n=1:numel(tt) t=tt(n); tp=t; shg subplot(121), cla draw_MBK2(xx(n,1),xx(n,2),u(t),Ms,Ks,Bs,w,'b'); axis([-1.5 1.5 -1 2]), axis('equal') subplot(122) plot(t,u(t),'r.', t,xx(n,1),'b.', t,xx(n,2),'g.') axis([0 tf -1 2]), hold on if n>1, pause(t-tp); end end end set(gcf,'color','white') shg function dx=dxdt(t,x) % State equation for a suspension system global M1 M2 K1 K2 B1 B2 %f=@(t)0; u=@(t)0.3*sin(4*t); du=@(t)1.2*cos(4*t); A=[ 0 0 1 0; 0 0 0 1; -K1/M1 K1/M1 -B1/M1 B1/M1; K1/M2 -(K1+K2)/M2 B1/M2 -(B1+B2)/M2]; dx=A*x+[0; 0; f(t)/M1; (K2*u(t)+B2*du(t)-f(t))/M2]; end function u=u(t) % Road disturbance u=0.2*sin(3*t); end function du=du(t) % Derivative of the road disturbance dt=1e-8; du=(u(t+dt)-u(t))./dt; end function f=f(t) f=10; end

function [ds,ls,Bs,Gammas]=imp_match_1stub_shunt(Z0,ZL,r,os) % Input: Z0 = Characteristic impedances [Z0,Z01] of the TL and stub % ZL = Load impedance % r = a+bi = Propagation constant % os = 'o'/'s' for open/short stub % Output: ds = Distance of the stub from load % ls = Length of the stub % Bs = Susceptance values to be cancelled out by the stub % Gammas = Intermediate points for ys=1+j*Bs if nargin<4, os='s'; end if imag(r)~=0, beta=imag(r); else beta=r; r=r*i; end lambda=2*pi/beta; if length(Z0)>1, Z01=Z0(2); Z0=Z0(1); else Z01=Z0; end yB=@(ZL,Z0,d)(Z0+ZL*tanh(r*d))./(ZL+Z0*tanh(r*d)); %Gamma=(1-yL)/(1+yL) % Voltage reflection coefficient %eq=@(d)real((Z0+ZL*tanh(r*d))/(ZL+Z0*tanh(r*d)))-1; % Eq.(7.3.5) eq=@(d)real(yB(ZL,Z0,d))-1; % Eq.(7.3.5) options=optimoptions('fsolve','Display','off','TolX',1e-20, 'TolFun',1e-20); ds=real([fsolve(eq,0,options) fsolve(eq,0.25*lambda,options)]); ds=mod(ds*beta,pi)/beta; %bs=imag((Z0+ZL*tanh(r*ds))./(ZL+Z0*tanh(r*ds)))*Z01/Z0; Bsr=Bs/Z0; bB=imag(yB(ZL,Z0,ds))*Z01/Z0; bs=-bB; Bs=bs/Z0; ys=1+j*bs; Gammas=(1-ys)./(1+ys); % Intermediate points %eq=@(d)imag((Z0+ZL*tanh(r*d))/(ZL+Z0*tanh(r*d)))+imag(yB); % such that Zi(d2)=1-j*gB %d2=real(fsolve(eq,0,options)) %ds=[d d2]; ds=ds+(ds<0)*pi/b; %bBs=[1 -1]*abs(imag(yB))*Z01/Z0; %yBs = 1+bBs*i; compare_with_Misra=1; if compare_with_Misra>0 YL=1/ZL; yL=YL*Z0; gL=real(yL); bL=imag(yL); A=gL*(gL-1)+bL^2; ds1=atan((bL+[-1 1]*sqrt(bL^2-A*(1-gL)))/A)/beta; % Eq.(5.1.2) of [M-1] ds1=mod(ds1*beta,pi)/beta; tbs=tan(beta*ds1); bs1=((bL+tbs).*(1-bL*tbs)-gL^2*tbs)./((gL*tbs).^2+(1-bL*tbs).^2)*Z01/Z0; % Eq.(5.1.3,4) of [M-1] [ds bs; ds1 bs1] end if lower(os(1))=='s', ls=acot(bB)/beta; % Eq.(7.3.7a) else ls=atan(-bB)/beta; % Eq.(7.3.7b) end % beta*(ls1+ls2)=pi/2 ls=mod(ls*beta,pi)/beta; % Usage: Z0=50; ZL=35-47.5i; lambda=1; r=2*pi*i/lambda; % [ds,ls]=imp_match_1stub_shunt(Z0,ZL,r,'s')

%pendulum2.m global G L1 M1 L2 M2 % Declare some global variables G=9.8; L1=2; M1=1; L2=1; M2=2; th1=120; w1=0; th2=-10; w2=0; % Set the initial values x0=deg2rad([th1; w1; th2; w2]); % Initial vector in rad tspan=[0 10]; % Time span of 10 sec % To solve the ODE [tt,xx] = ode45(@dxdt,tspan,x0); x1 = L1*sin(xx(:,1)); y1 = -L1*cos(xx(:,1)); x2 = L2*sin(xx(:,3))+x1; y2 = -L2*cos(xx(:,3))+y1; % Animation figure(1), calf plot(0,0,'o'); hold on ax=gca; % Current figure axis axis('equal') ax.XLim=[-3.5 3.5]; ax.YLim=[-3.5 3.5]; for n=1:numel(tt) cla % Clear the current figure plot(0,0,'o','MarkerSize',2); plot([0 x1(n)],[0 y1(n)],'r') % Rod 1 plot(x1(n),y1(n),'ro','MarkerSize',M1*5) plot([x1(n) x2(n)],[y1(n) y2(n)],'g') % Rod 2 plot(x2(n),y2(n),'go','MarkerSize',M2*5) text(-0.6,2.5,"Timer: "+num2str(tt(n),2)); if n>1, pause(tt(n)-tt(n-1)); end end set(gcf,'color','white') function dx=dxdt(t,x) % State equation for a double pendulum global G L1 M1 L2 M2 th1=x(1); w1=x(2); th2=x(3); w2=x(4); s1=sin(th1); s2=sin(th2); s12=sin(th1-th2); c12=cos(th1-th2); A = [L1*(M1+M2) L2*M2*c12; L1*M2*c12 L2*M2]; %Eq. (P9.14.3) b = [-L2*M2*s12*w2^2-G*(M1+M2)*s1; L1*M2*s12*w1^2-G*M2*s2]; dx24 = A\b; dx = [w1; dx24(1); w2; dx24(2)]; end % Reference 1: "Engineering Mathematics with MATLAB" by Won Y. Yang et al. (www.google.co.kr/books/edition/_/R2OuEAAAQBAJ?hl=ko&gbpv=1&pg=PR1&dq=Engineering+Mathematics+with+MATLAB) % Reference 2: "Applied Numerical Methods Using MATLAB" by Won Y. Yang et al. (www.google.co.kr/books/edition/Applied_Numerical_Methods_Using_MATLAB/LcnYDwAAQBAJ?hl=ko&gbpv=1&dq=Applied+Numerical+Methods+Using+MATLAB&printsec=frontcover)

amazing! it can be better with a voice explaining!

Very nice 😮

교수님! 은퇴하신 후에도 후학들과 교육환경 개선을 위해 노력하시는 열정에 경의를 표합니다. 교수님의 선한 생각이 널리퍼져서 더 나은 대한민국이 되길 희망합니다

%Here is the MATLAB function saved in the m-file named "pde_wave.m", that can be downloaded at <onedrive.live.com/?authkey=%21AE%2DoDXK83bmhRVY&id=85B008EBE473BF17%2130172&cid=85B008EBE473BF17>. function [u,x,t]=pde_wave(a,xf,tf,it0,i1t0,bx0,bxf,M,N) %solve a^2*u_xx = u_tt = for 0<=x<=xf, 0<=t<=tf % Initial Condition: u(x,0) = it0(x), u_t(x,0) = i1t0(x) % Boundary Condition: u(0,t) = bx0(t), u(xf,t)=bxf(t) % M = Number of subintervals along x-axis % N = Number of subintervals along t-axis dx= xf/M; x= [0:M]'*dx; dt= tf/N; t= [0:N]*dt; for i=1:M+1, u(i,1)= it0(x(i)); end for k=1:N+1 u([1 M+1],k)= [bx0(t(k)); bxf(t(k))]; end r= (a*dt/dx)^2, r1= r/2; r2= 2*(1-r); u(2:M,2)= r1*u(1:M-1,1) +(1-r)*u(2:M,1) +r1*u(3:M+1,1) +dt*i1t0(x(2:M)); for k=3:N+1 u(2:M,k)= r*u(1:M-1,k-1) +r2*u(2:M,k-1) +r*u(3:M+1,k-1) -u(2:M,k-2); end