גרפים: הבדלים בין גרסאות
מתוך מעבדת מבוא בחשמל
(הוספת הערה לשאלות מסוג מאקסימה) |
(הורדת קוד מיותר בשליפת תדר בסריקה ליניארית) |
||
שורה 1,090: | שורה 1,090: | ||
/* calculate f according to time */ | /* calculate f according to time */ | ||
f: (fT-f0)*t/T+f0; | f: (fT-f0)*t/T+f0; | ||
− | |||
− | |||
− | |||
/* do the plot */ | /* do the plot */ | ||
− | p: plot(v_out, [t, 0, T], [ylabel, "Voltage [V]"], [xlabel, "Timing [S]"], grid2d, [y, 0, volt_size], [xtics, tC, | + | p: plot(v_out, [t, 0, T], [ylabel, "Voltage [V]"], [xlabel, "Timing [S]"], grid2d, [y, 0, volt_size], [xtics, tC, T, T], [ytics, 0, v_cutoff, v_cutoff], [label, ["0", 0, 0], [string(T), T, 0]]); |
</syntaxhighlight> | </syntaxhighlight> | ||
גרסה מתאריך 21:28, 26 בדצמבר 2018
תוכן עניינים
- 1 מטלאב
- 1.1 גרף הדגמת FFT ו-THD
- 1.2 ייצוג ספקטראלי של גלים יסודיים
- 1.3 בניית גלים יסודיים באמצעות הרמוניות
- 1.4 גרף הפרש-מופע
- 1.5 גרף עקומות ליסאז'ו
- 1.6 תגובה לתדר של מעגל RC
- 1.7 אופיין חשל מגנטי
- 1.8 גרף תופעת מעבר במעגל RL
- 1.9 תגובה לתדר ותדר-זוויתי של מעגל RL טורי
- 1.10 תופעת מעבר במעגל RC עם מקור סינוסי
- 1.11 גרף מצבי ריסון
- 1.12 ניתוח נתונים מאורקד
- 1.13 שליפת הפרש-מופע מעקומת ליסאז'ו
- 2 ג'אווה-סקריפט
- 3 מאקסימה
1 מטלאב
1.1 גרף הדגמת FFT ו-THD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FFT with THD example %%%
%%% based on https://www.mathworks.com/help/matlab/ref/fft.html %%%
%%% with a few modifications %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 300; % Length of signal
t = (0:L-1)*T; % Time vector
freq = 50; % frequency
S1 = 5*sin(2*pi*freq*t); % 1st (base) harmony
S2 = 2*sin(2*pi*3*freq*t); % 2nd harmony distortion
S = S1 + S2;
f = Fs*(0:(L/2))/L;
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
figure('Name','THD example');
subplot(311);
plot(t,S1,t,S2);
title('Two pure sine waves');
ylabel('Amplitude [V]');
xlabel('Time [S]');
legend(['1st harmony: ',num2str(freq),' Hz'],['2nd harmony: ',num2str(3*freq),' Hz']);
subplot(312);
plot(t,S);
title(['Addition of the above two sine waves with frequencies of ',num2str(freq),' and ',num2str(3*freq),' Hz']);
ylabel('Amplitude [V]');
xlabel('Time [S]');
subplot(313);
plot(f,P1);
title('Single-Sided Amplitude Spectrum of 1st (base) harmony + 2nd harmony');
xlabel('f (Hz)');
ylabel('Amplitude [V]');
set(gcf, 'Position', [600, 80, 800, 900]);
1.2 ייצוג ספקטראלי של גלים יסודיים
יש לשנות את המשתנה wave לערכים של בין 1 ל-6 כדי לקבל את ששת הגלים השונים:
- סינוסי
- ריבועי
- שן-מסור
- משולש
- חצי-גל סינוסי מיושר
- סדרת-דפקים
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FFT based on %%%
%%% https://www.mathworks.com/help/matlab/ref/fft.html %%%
%%% with a few modifications %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
freq = 1000; % frequency
Fs = freq*20; % Sampling frequency
T = 1/Fs; % Sampling period
L = 200; % Length of signal
t = (0:L-1)*T; % Time vector
% three different waves
RMS_Amplitude = 1;
x = 2*pi*freq*t;
d=20/100; % Duty cycle for pulse wave
Names = {'Sine', 'Square', 'Sawtooth', 'Triangle','Half rectified sine','Pulse'};
Amplitudes = [sqrt(2), 1, sqrt(3), sqrt(3), 2, 1/sqrt(d)] * RMS_Amplitude;
half_wave = sin(x); half_wave(find(half_wave<0))=0;
pulse_wave = pulstran(x,(0:9)*max(x)/10,@rectpuls,d*max(x)/10);
Waves = {sin(x), square(x), sawtooth(x), sawtooth(x,0.5), half_wave, pulse_wave};
wave = 6; % actual wave to display
S=Amplitudes(wave)*Waves{wave};
f = Fs*(0:(L/2))/L;
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
figure('Name','FFT example');
set(gcf, 'Position', [600, 80, 800, 900]);
ax=subplot(211);
plot(t,S); % plot wave
title(['Original ',Names{wave},' wave']);
ylabel('Amplitude [V]');
xlabel('Time [S]');
if max(ax.YLim)==max(S)
ax.YLim=ax.YLim*1.1;
end
ax=subplot(212);
P1_RMS=P1;
non_DC=find(f>freq/2); % find non DC frequencies
P1_RMS(non_DC)=P1_RMS(non_DC)/sqrt(2); % calculate RMS
plot(f,P1_RMS); % plot FFT of wave
title(['FFT of the ',Names{wave},' wave']);
xlabel('f (Hz)');
ylabel('Amplitude [RMS V]');
if max(ax.YLim)==max(P1_RMS)
ax.YLim=ax.YLim*1.1;
end
1.3 בניית גלים יסודיים באמצעות הרמוניות
זוהי תוכנה הופכית ממה שמופיע בייצוג ספקטראלי של גלים יסודיים, הפעם באמצעות נוסחא של ההרמוניות - בונים את הגל הסופי.
יש לשנות את המשתנה wave_num לערכים של בין 1 ל-6 כדי לקבל את ששת הגלים השונים:
- גל ריבועי
- סדרת-דפקים
- מיישר חצי-גל
- מיישר גל-מלא
- גל שן-מסור
- גל משולש
% structure for fourer waves
% 1: Name, 2: RMS Coefficient, 3: DC Average, 4: Series of harmonies, 5: Stand alone harmony if any
fourier_waves = {;
{'Square', '1', '0', '4*A/pi*sin((2*n-1)*w0.*t)/(2*n-1)'};
{'Pulse', '1/sqrt(d)', 'A*d', '2*A/(n*pi)*sin(n*pi*d)*cos(n*w0.*t)'};
{'Half wave rectified sine', '2', 'A/pi', '-2*A/pi*cos(2*n*w0.*t)/(4*n^2-1)', 'A/2*sin(w0.*t)'};
{'Full wave rectified sine wave', 'sqrt(2)', '2*A/pi', '-4*A/pi*cos(2*n*w0.*t)/(4*n^2-1)'};
{'Saw tooth wave', 'sqrt(3)', '0', '-2*A/pi*sin(n*w0.*t)/n'};
{'Triangle wave', 'sqrt(3)', '0', '8*A/pi^2*cos((2*n-1)*w0.*t)/(2*n-1)^2'}
};
wave_num = 2; % which wave to depict
wave = fourier_waves{wave_num}; % get data for selected wave
A_RMS = 1; % Amplitude in volts RMS
d=20/100; % Duty cycle for pulse wave
RMS_coef=eval(wave{2});
A=A_RMS * RMS_coef;
f = 1000; % frequency in Hz
harmonies = 1000; % How many harmonies to sum up
T=1/f; % period
t=linspace(-T,T,1001); % resolution of two full periods
w0=2*pi*f; % fundamental harmony
y=zeros(size(t)); % allocate memory for the harmonics addition
figure('Name', [wave{1},' wave']); % create a plot window
set(gcf, 'Position', [600, 80, 800, 900]); % change window size
fourier_wave_prefix=wave{3}; % prefix of Fourier series
fourier_wave=wave{4}; % addends of the Fourier series
if length(wave)==5
fourier_wave1=wave{5}; % out of serise harmony
else
fourier_wave1=0;
end
subplots=9; % how many subplot steps to show
rows=sqrt(subplots); cols=rows;
axs=cell(1,subplots); % save all axes for normalization
y_max=0; y_min=0;
for harmony=0:harmonies % sums up all the harmonis
switch harmony
case 0
n=0;
% calculate DC (harmony 0)
f_n=eval(fourier_wave_prefix);
case 1
% calculate another harmony
if fourier_wave1==0
f_n=eval(fourier_wave);
else
f_n=eval(fourier_wave1);
n=n-1;
fourier_wave1=0;
end
otherwise
f_n=eval(fourier_wave);
end
y=y+f_n;
% plot the intermediary series summation
if harmony<subplots-1
axs{harmony+1}=subplot(rows,cols,harmony+1);
plot(t,y);
switch harmony
case 0
title('DC (harmony 0)');
case 1
title('First harmony only');
otherwise
title(sprintf('Summation of %d harmonies',harmony));
end
y_min=min(min(y), y_min); y_max=max(max(y), y_max);
end
n=n+1;
end
axs{subplots}=subplot(rows,cols,subplots);
plot(t,y); % plot the final wave
title(['Summation of ',num2str(harmonies),' harmonies']);
% normalize all axes
dy=y_max-y_min;
ylim=[y_min-dy/10, y_max+dy/10];
for i=1:subplots
axs{i}.YLim=ylim;
end
% plots actual wave (with 1k harmonies) in a new window
figure('Name','Actual wave');
plot(t,y);
ax=gca;
ax.YLim=ylim;
xlabel('Time [S]');
ylabel('Amplitude [V]');
title([wave{1},' wave']);
1.4 גרף הפרש-מופע
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plots two sine waves with %%%
%%% phase and amplitude differences %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = figure; % create a plot window
cycles = 2+1/16; % display about 2 cycles
frequency = 1000; % 1 kHz frequency
total_x = 2*pi*cycles; % total horizontal length
x = linspace(0,total_x,1000); % horizontal axis - time scale
A1 = 1; % 1st amplitude
y1 = A1*sin(x); % vertical axis - voltage/current scale
phase = 2*pi/3; % create a phase difference
A2 = 0.7; % 2nd amplitude
y2 = A2*sin(x-phase); % another wave with phase difference
plot(x,y1,'k',x,y2,'b','LineWidth',3); % plot the graphs
h1 = gca; % Get handle to Current Axis (1)
h1.YLim = [-1.25, 1.25]; % set scale of figure 1
left1 = pi/2; % left cursor at maximum point
right1 = left1 + 2*pi; % next peak
dx = (right1-left1)/20;
annotation('doublearrow',x_to_norm_v2(left1+dx, right1-dx),y_to_norm_v2(A1, A1),'Color','k','LineWidth',2);
text(left1+(right1-left1)/2,A1+0.1,'T','Color','k','HorizontalAlignment','Center');
refline(0,0); % add horizontal reference line
left2 = pi/2 + phase;
right2 = left2 + 2*pi; % next peak
annotation('doublearrow',x_to_norm_v2(left2+dx, right2-dx),y_to_norm_v2(A2, A2),'Color','b','LineWidth',2);
text(left2+(right2-left2)/2,A2+0.1,'T','Color','b','HorizontalAlignment','Center');
% plot vertical lines (cursors) for peaks
hold on;
cursor1x=[left1, left1]; cursor1y=[-0.75, A1]; % left cursor
cursor2x=[left2, left2]; cursor2y=[-0.75, A2]; % right cursor
plot(cursor1x, cursor1y, 'k--', cursor2x, cursor2y, 'b--','LineWidth',2); % plot cursors
annotation('doublearrow',x_to_norm_v2(left1+dx, left2-dx), y_to_norm_v2(-0.75, -0.75), 'Color','r','LineWidth',2); % add ruler
text(left1+(left2-left1)/2,-0.75+0.1,'\Delta T','Color','r','HorizontalAlignment','Center');
% plot vertical lines (cursors) for zero crossings
hold on;
left1z=2*pi;
right1z=2*pi+phase;
cursor1xz=[left1z, left1z]; cursor1yz=[-1, 0]; % left cursor
cursor2xz=[right1z, right1z]; cursor2yz=[-1, 0]; % right cursor
plot(cursor1xz, cursor1yz, 'k--', cursor2xz, cursor2yz, 'b--','LineWidth',2); % plot cursors
annotation('doublearrow',x_to_norm_v2(left1z+dx, right1z-dx), y_to_norm_v2(-1, -1), 'Color','m','LineWidth',2); % add ruler
text(left1z+(right1z-left1z)/2,-1+0.1,'\Delta T','Color','m','HorizontalAlignment','Center');
1.5 גרף עקומות ליסאז'ו
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plots Lissajous curve shapes %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('Name','Liss');
% draw a unit circle in the middle
subplot(5,5,[7:9,12:14,17:19]);
phi_ind=0:11;
phi_rad=phi_ind*pi/6; % phase shifts in radians
phi_deg=phi_ind*30; % phase shifts in degrees
z=exp(j*phi_rad);
compass(z);
numbers=findall(gcf,'Type','text');
for i=1:numel(numbers)
txt = numbers(i).String;
if isempty(regexp(txt,'0$'))
% erase non degrees numbers
numbers(i).Visible='off';
else
% add degrees symbol
numbers(i).String=[txt,'^\circ'];
end
end
xt=linspace(0,2*pi,1000);
x=sin(xt);
plot2angle=[15, 10, 4, 3, 2, 6, 11, 16, 22, 23, 24, 20];
for i=phi_ind
ax = subplot(5,5,plot2angle(i+1));
y=sin(xt + phi_rad(i+1)); % add phase to y
plot(x,y); % plot the curve
% increase gap between plot and frame
ax.XLim=ax.XLim*1.1;
ax.YLim=ax.YLim*1.1;
% remove labels
ax.YTickLabel={};
ax.XTickLabel={};
% create title for each subplot
phis = [phi_deg(i+1), phi_deg(i+1)-360];
if (phis(1) > 180)
phis = fliplr(phis);
else
if phis(1) == 0
phis(2) = 360; % instead of -360
end
end
title(sprintf('%d^\\circ = %d^\\circ', phis));
grid on;
end
% equalize the scale of y and x axes
set(gcf, 'Position', [100, 100, 900, 1000]);
1.6 תגובה לתדר של מעגל RC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Create a frequency response for RC circuit %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = 1000; % 1kHz frequency
R = 1000; % 1kΩ Resistance
C = 1/(2*pi*f*R); % 0.159µF Capacitance
numer = [1]; denom = [R*C*1 1]; % Hc(s)=1/(1+RCs)
w = logspace(0,6,1000); % 10^0 to 10^6 angular-frequency range
freqs(numer,denom,w); % create the graph
set(gcf, 'Position', [600, 80, 800, 900]); % enhance resolution
1.7 אופיין חשל מגנטי
%%%%%%%%%%%%%%%%%%%%%%
%%% Plot B-H curve %%%
%%%%%%%%%%%%%%%%%%%%%%
try
% worked on matlab version R2017a
bhsim = sim('elec_inductor_hysteresis','MaxStep','1e-5');
simlog = get(bhsim,'simlog_elec_inductor_hysteresis');
catch
% worked on matlab version R2016b
bhsim = sim('elec_inductor_with_hysteresis','MaxStep','1e-5');
simlog = get(bhsim,'simlog');
end
inductor=simlog.Nonlinear_Inductor.inductor;
H=inductor.H.series.values;
B=inductor.B.series.values;
figure('Name','B-H Hysteresis');
plot(H,B);
set(gcf, 'Position', [100, 100, 900, 1000]);
grid on;
xlabel('H [A/m]');
ylabel('B [T]');
% to disable 'Error due to multiple causes' for next runs
close_system();
1.8 גרף תופעת מעבר במעגל RL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Draw transient reponse of RL circuit %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V=10; % power-supply voltage [V]
R=10000; % resistor value [Ω]
L=10; % inductor value [H]
tau=L/R; % time constant tau [s]
t=linspace(0,10*tau,1001); % time axis
il=V/R+(0-V/R)*exp(-t/tau);
figure('Name','RL circuit');
plot(t,il,'LineWidth',2);
% signage
xlabel('time [s]');
ylabel('i_L [A]');
title('RL circuit transient response');
hold on;
% extract special tau points
x = (1:5)*tau;
[~,il_tau]=intersect(t,x);
y = il(il_tau);
% calculate axis limits
max_y = max(y) * 1.2;
max_x = max(t) * 1.1;
min_y = -max_y / 12;
min_x = -max_x / 11;
% draw cross lines to emphasize special tau points
for i = 1:length(x)
plot([min_x,x(i),x(i)],[y(i),y(i),min_y]);
end
legend('i_L','1\tau','2\tau','3\tau','4\tau','5\tau');
% change axis
ax=gca;
ax.YLim=[min_y max_y];
ax.YAxis.Exponent = -3;
ax.XLim=[min_x max_x];
1.9 תגובה לתדר ותדר-זוויתי של מעגל RL טורי
R=10000; % Resistor value [Ω]
L=10; % Inductor value [H]
tau=L/R; % time constant
cutoff=1/(2*pi*tau); % Cutoff frequency [Hz]
f=linspace(0,cutoff*5,1000); % Frequency sweep [Hz]
w=2*pi*f; % Angular frequency sweep [rad/sec]
Vin=10;
ZL=j*w*L; % Inductor Impedance
Vr=Vin*R./(R+ZL); % Complex Voltage of resistor
VL=Vin*ZL./(R+ZL); % Complex Voltage of inductor
mag_vr=abs(Vr);% get the magnitude of resistor's voltage
mag_vL=abs(VL);% get the magnitude of inductor's voltage
% Plot voltage vs. angular frequency
figure; hold on;
set(0, 'DefaultLineLineWidth', 2);
plot(w,mag_vr,w,mag_vL); % plot graphs
plot([2*pi*cutoff 2*pi*cutoff],[0 Vin]); % mark the cutoff vertically
plot([0 max(w)],[Vin/sqrt(2) Vin/sqrt(2)]); % mark the cutoff horizontally
% add signage
title('V_{OUT} vs. Angular Frequency');
xlabel('Angular Frequency [rad/s]');
ylabel('Voltage [V]');
legend('Resistor voltage','Inductor voltage','Cutoff X','Cutoff Y');
% Plot voltage vs. frequency
figure; hold on;
plot(f,mag_vr,f,mag_vL); % plot graphs
plot([cutoff cutoff],[0 Vin]); % mark the cutoff vertically
plot([0 max(f)],[Vin/sqrt(2) Vin/sqrt(2)]); % mark the cutoff horizontally
% add signage
title('V_{OUT} vs. Frequency');
xlabel('Frequency [1/s]=[Hz]');
ylabel('Voltage [V]');
legend('Resistor voltage','Inductor voltage','Cutoff X','Cutoff Y');
1.10 תופעת מעבר במעגל RC עם מקור סינוסי
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plot RC circuit transient response with Sine wave input %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta=pi/4; % Phase angle of AC input
AC = calcRC(theta,'ac');
% RC Transient response equation
RCTrans=calcRC(theta,'trans');
RCSteadyState=calcRC(theta,'steady');
RCTransWithSine=RCTrans+RCSteadyState;
% calculate exponential decay
RCDecay=calcRC(theta,'decay');
% calculate theta which will cancel the transient response
corrected_theta=calcRC(0,'correction');
RCNoTrans=calcRC(corrected_theta,'both');
% plot input and output of original theta
figure('Name','RC with AC source');
subplot(2,1,1);
t=calcRC(0,'time');
plot(t,AC,t,RCTransWithSine,t,RCDecay,[min(t),max(t)],[0,0]);
legend('AC Input','Output (capacitor)','Decay envelope','Reference point');
xlabel('Time [s]'); ylabel('Voltage [v]');
title(calcRC(theta,'title'));
% plot input and output of corrected theta
subplot(2,1,2);
plot(t,calcRC(corrected_theta,'ac'),t,RCNoTrans,[min(t),max(t)],[0,0]);
legend('AC input shifted','Output (capacitor)','Reference point');
xlabel('Time [s]'); ylabel('Voltage [v]');
title(calcRC(corrected_theta,'title'));
% change axis
ax=gca;
top=max(AC);
ax.YLim=[-top-1, top+1];
set(gcf, 'Position', [100, 100, 600, 800]);
% AC input/Transient/SteadyState response of RC circuit
function response = calcRC(angle,role)
Vm=10; % Amplitude of sine wave
V0=2; % Initial voltage of capacitor
f=700; % frequency of AC input
w=2*pi*f; % angular frequency
count=10; % how many time constants to plot
R=10e3; % Resistor of 10kΩ
C=0.1e-6; % Capacitor of 0.1µF
tau=R*C; % Time constant
t=linspace(0,count*tau,1000); % Time axis
switch lower(role)
case 'title'
response = sprintf('V_S(t)=%dsin(2\\pi\\cdot%dt+%d^\\circ); V_0=%dV',Vm,f,round(angle/pi*180),V0);
case 'time'
%%% returns the time axis
response = t;
case 'correction'
%%% returns angle to cancel out transient response
response = asin(V0/Vm*sqrt(1+(w*tau)^2))-atan(-w*tau);
case 'both'
%%% return entire response (transient and steady state)
response = calcRC(angle,'trans') + calcRC(angle,'steady');
case 'trans'
%%% return transient response
response = (V0-Vm*sin(angle+atan(-w*tau))/sqrt(1+(w*tau)^2))*exp(-t/tau);
case 'steady'
%%% return steady state response
response = Vm/sqrt(1+(w*tau)^2)*sin(w*t+angle+atan(-w*tau));
case 'ac'
%%% AC input of circuit;
response = Vm*sin(w*t+angle);
case 'decay'
%%% returns decay envelope
% max voltage of full response
[peakVolt, peakInd]=max(calcRC(angle,'both'));
peakTime=t(peakInd);
% voltage diff of peak and steady state
MAXSteady=max(calcRC(angle,'steady'));
dV=peakVolt-MAXSteady;
response = dV*exp(-(t-peakTime)/tau)+MAXSteady;
end
end
1.11 גרף מצבי ריסון
%%%%%%%%%%%%%%%%%%%%%%
%%% Damping ratios %%%
%%%%%%%%%%%%%%%%%%%%%%
figure('Name','Damping');
hold on;
% iterate thru a few different resistor values
Rs=[1000,5000,20e3,100e3,1e6];
Rs_len = length(Rs);
h = zeros(1,Rs_len);
legend_text = cell(1,Rs_len);
t=linspace(0,0.005,10000);
for i=1:Rs_len
R=Rs(i); % Current resitor in Ω
vc=solve_ode(R); % solve current ODE
vct=double(vc(t)); % fill in values
h(i)=plot(t,vct,'LineWidth',2); % plot current solution
% create description for each resistor type
if R < 5e3
info = '- overdamped';
elseif R == 5e3
info = '- critically damped';
elseif R >= 1e6
info = '~ undamped (lossless)';
else
info = '- underdamped';
end
legend_text{i}=['R=',num2str(R),' \Omega ',info];
legend(h(1:i),legend_text{1:i}); % update legend
drawnow; % view plots during run
end
% signage
xlabel('Time [S]');
ylabel('V_C [V]');
title('Different damping ratios');
grid on;
set(gcf,'Position',[100,100,800,600])
% solve second order ode of type
% y''+2α*y'+w^2*y=w^2*V
% y(0)=V0
% y'(0)=-1/C*(V0/R+I0);
function solution = solve_ode(R)
C=0.01e-6; % Capacitor in F
L=1; % Inductor in H
tau=R*C; % Time constant
alpha=1/(2*tau); % Damping factor
w0=1/sqrt(L*C); % natural frequency
V = 10; % Input voltage
V0=5; % Initial capacitor voltage
I0=0; % Initial inductor current
dVc0=-1/C*(V0/R+I0); % Derivate of capacitor's voltage at t=0
% Based on
% https://www.mathworks.com/help/symbolic/solve-a-single-differential-equation.html
syms y(t)
Dy = diff(y);
ode = diff(y,t,2)+2*alpha*Dy+w0^2*y == w0^2*V;
cond1 = y(0) == V0;
cond2 = Dy(0) == dVc0;
conds = [cond1, cond2];
ySol(t) = dsolve(ode, conds);
solution = simplify(ySol);
end
1.12 ניתוח נתונים מאורקד
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analyze rectifier data from orcad %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
csv = load('orcaddata.csv');
% get timestamps
time = csv(:,1);
% get voltage of resistor
vr = csv(:,2);
% maximum voltage
Vm = max(vr);
% calculate average
VDC = mean(vr);
% calculate rms
VDCRMS = rms(vr);
% calculate form factor
FF = VDCRMS/VDC;
% calculate ripple factor
RF = sqrt(FF^2-1);
% plot all data
figure('Name','Rectifier');
total_x=[0 max(time)];
plot(time,vr,total_x,[Vm Vm],total_x,[VDC VDC],0,0,0,0,0,0);
txtVm=['V_m=',num2str(Vm),'V'];
txtVDC=['V_{DC}=',num2str(VDC),'V'];
txtVDCRMS=['V_{DCRMS}=',num2str(VDCRMS),'V'];
txtFF=['Form Factor=',num2str(FF)];
txtRF=['Ripple Factor=',num2str(RF)];
legend('V_R',txtVm,txtVDC,txtVDCRMS,txtFF,txtRF);
xlabel('Time [s]');
ylabel('Amplitude [V]');
ax=gca;
ax.YLim=[min(vr)-1, max(vr)+1];
1.13 שליפת הפרש-מופע מעקומת ליסאז'ו
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plot of Lissajous curve with %%%
%%% markes to extract the phase %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = linspace(0,2*pi,1000); % time between 0 and 2*pi
f = 1; % frequency in Hz
V1 = 4; % amplitude of first wave
wave1 = V1*sin(2*pi*f*t); % a sine wave
phi = pi/4; % phase difference of 45 degrees
V2 = 30; % amplitude of second wave
wave2 = V2*sin(2*pi*f*t+phi); % same sine wave but with a phase
figure('Name','Lissajous');
hold on;
plot(wave1, wave2, 'k','LineWidth',3); % plot the graphs
ax = gca; % Get handle of Current Axis (1)
ax.YLim = ax.YLim*1.25; % increase Y scale
ax.XLim = ax.XLim*1.25; % increase X scale
% calculate the contact points on the Lissajous curve
a1 = [0, abs(V2*sin(phi))]; % top a
a2 = [0, -abs(V2*sin(phi))]; % bottom a
b1 = [V1*cos(phi), abs(V2)]; % top b
b2 = [-V1*cos(phi), -abs(V2)]; % bottom b
% write the names of the points
gap_x = V1*0.1;
gap_y = V2*0.1;
text(a1(1)+gap_x, a1(2), 'a1', 'FontSize', 14, 'Color', 'r');
text(a2(1)+gap_x, a2(2), 'a2', 'FontSize', 14, 'Color', 'r');
text(b1(1), b1(2)+gap_y, 'b1', 'FontSize', 14, 'Color', 'b');
text(b2(1), b2(2)-gap_y, 'b2', 'FontSize', 14, 'Color', 'b');
% plot horizontal lines for the contact points
a1x=[-V1/2, a1(1)]; a1y=[a1(2), a1(2)]; % top middle
plot(a1x, a1y, 'r--', 'LineWidth',2); % plot line
a2x=[-V1/2, a2(1)]; a2y=[a2(2), a2(2)]; % bottom middle
plot(a2x, a2y, 'r--', 'LineWidth',2); % plot line
b1x=[0, b1(1)+gap_x]; b1y=[b1(2), b1(2)]; % top middle
plot(b1x, b1y, 'b--', 'LineWidth',2); % plot line
b2x=[b2(1), b1(1)+gap_x]; b2y=[b2(2), b2(2)]; % bottom middle
plot(b2x, b2y, 'b--', 'LineWidth',2); % plot line
% plot a vertical line between a1 and a2
annotation('doublearrow',x_to_norm_v2(-V1/3, -V1/3),y_to_norm_v2(a1(2), a2(2)),'Color','r','LineWidth',2);
% plot a vertical line between b1 and b2
annotation('doublearrow',x_to_norm_v2(b1(1), b1(1)),y_to_norm_v2(b1(2), b2(2)),'Color','b','LineWidth',2);
% write the names of the vertical lines
text(-V1/3+gap_x, a1(2)/2, 'a', 'FontSize', 14, 'Color', 'r');
text(b1(1)+gap_x, -V2/2, 'b', 'FontSize', 14, 'Color', 'b');
% turn on grid and draw X, Y axes
grid on;
gray=0.4*ones(1,3);
annotation('arrow',x_to_norm_v2(ax.XLim(1), ax.XLim(2)),y_to_norm_v2(0, 0),'Color',gray,'LineWidth',2);
annotation('arrow',x_to_norm_v2(0, 0),y_to_norm_v2(ax.YLim(1), ax.YLim(2)),'Color',gray,'LineWidth',2);
text(0-gap_x,V2+gap_y,'y', 'FontSize', 14, 'Color', gray);
text(V1+gap_x,0-gap_y,'x', 'FontSize', 14, 'Color', gray);
% draw green circles as the points
for dots={a1,a2,b1,b2}
dot=dots{1};
plot([dot(1), dot(1)], [dot(2), dot(2)], 'o', 'Color', 'g', 'MarkerSize',10, 'MarkerFaceColor','g');
end
2 ג'אווה-סקריפט
2.1 מסנן מעביר גבוהים או נמוכים
var rand = {random};
// high pass = 1, low pass = 0
var high_pass = Math.round(rand);
// grid frequency (MUST BE ROUND)
var grid_frequency = Math.round(6*rand+4);
// maximum voltage
var volt_size = Math.round((10*rand+2)*100)/100;
// Maximum frequency for sweep is between 100 kHz and 10 MHz
var max_frequency = Math.pow(10,Math.round(3*rand+5));
// frequency of internal sine wave
var frequency = 30;
// the frequency in the middle of the sweep
var middle_frequency = max_frequency / 2;
// Resistor size in Ohm
var R = 160;
// Capcitor size in Farad
var C = 1e-6;
// define a canvas dimensions along with grid lines
var w=500;
var h=360;
var grid_y=grid_frequency;
var grid_x=grid_frequency*2-1;
var hs=h*(grid_y-1)/grid_y; // sub height
// create a canvas inside the table
var c = document.createElement("canvas");
c.setAttribute("width",w);
c.setAttribute("Height",h);
var element = document.getElementById("insertcanvashere_sweep");
element.appendChild(c);
var ctx = c.getContext("2d");
// set style of text
ctx.font="20px Georgia";
ctx.save();
// set the style of grid lines
ctx.setLineDash([1, 3]);/*dashes are 5px and spaces are 3px*/
ctx.strokeStyle="blue";
// dispersion of x grid lines
var skip_x = w/(grid_frequency*2);
var y;
// draw x grid
var half = Math.floor(grid_y/2);
// pixels amount of each y gridline
var skip_h = h/grid_frequency;
// position in which to right voltage
var mark_y = Math.round(rand*(half-2))+1;
for (var i=0; i<grid_y; i++) {
y=(-half+i)*skip_h+h/2;
ctx.beginPath();
ctx.moveTo(0,y);
ctx.lineTo(w,y);
ctx.stroke();
// write voltages
if (i==half) {
ctx.fillText("0 V",w/2,y - 5);
// draw different color for the period
ctx.save();
ctx.setLineDash([1, 0]);
ctx.strokeStyle="black";
ctx.stroke();
ctx.restore();
}
if (i==mark_y) {
// calculate voltage at current point
var voltage = (half - mark_y) * skip_h * 2 / hs * volt_size;
round_volt = Math.round(voltage * 1000) / 1000;
ctx.fillText(round_volt.toString() + 'V',w/2,y - 5);
// draw different color for the period
ctx.save();
ctx.setLineDash([10, 6]);
ctx.strokeStyle="black";
ctx.stroke();
ctx.restore();
}
}
min_y=y;
// draw y grid
var period=middle_frequency / grid_frequency;
var total_period=0;
var x_counter=0;
for (var x=skip_x; x<w; x+=skip_x) {
ctx.beginPath();
ctx.moveTo(x,0);
ctx.lineTo(x,h);
ctx.stroke();
// write time period
total_period+=period;
x_counter++;
if ((x_counter==2) || (x_counter==grid_x-2)) {
var round_period = Math.round(total_period * 10) / 10;
var units = 'Hz';
if (x_counter != 2) {
round_period = total_period / 1000;
round_period = Math.round(round_period * 10) / 10;
units = 'kHz';
}
ctx.fillText(round_period.toString() + units,x,h - 10);
// draw different color for the period
ctx.save();
ctx.setLineDash([10, 6]);
ctx.strokeStyle="black";
ctx.stroke();
ctx.restore();
}
}
ctx.restore();
// choose LPF
var filter_function='1/(omega*R*C+1)';
if (high_pass == 1) {
// choose HPF
filter_function='omega*R*C/(omega*R*C+1)';
}
// draw a positive RC response graph
ctx.beginPath();
ctx.moveTo(0,h/2);
var max_frequency_log=Math.log10(max_frequency);
for (var x=1; x<w; x++) {
// calculate the logarithmic frequency
current_frequency=Math.pow(10,x/w*max_frequency_log);
var omega = 2 * Math.PI * current_frequency;
y=eval(filter_function);
y=h/2-y*hs/2;
ctx.lineTo(x,y);
}
ctx.strokeStyle="red";
ctx.stroke();
// draw a negative RC response graph
ctx.beginPath();
ctx.moveTo(0,h/2);
for (var x=1; x<w; x++) {
// calculate the logarithmic frequency
current_frequency=Math.pow(10,x/w*max_frequency_log);
var omega = 2 * Math.PI * current_frequency;
y=eval(filter_function);
y=y*hs/2+h/2;
ctx.lineTo(x,y);
}
ctx.strokeStyle="red";
ctx.stroke();
// draw a sine wave inside the RC response
ctx.beginPath();
ctx.moveTo(0,h/2);
for (var x=1; x<w; x++) {
current_frequency=Math.pow(10,x/w*max_frequency_log);
var omega = 2 * Math.PI * current_frequency;
var amplitude=eval(filter_function);
y=amplitude*hs/2*Math.sin(2*Math.PI*frequency*x/w)+h/2;
ctx.lineTo(x,y);
}
ctx.strokeStyle="green";
ctx.stroke();
2.2 ניתוח גל סינוס - הוצאת תדר
var total_time = {total_time};
var volt_size = {volt_size};
var frequency = {frequency};
var w=500;
var h=360;
var grid_y=7;
var grid_x=frequency*2-1;
var hs=h*(grid_y-1)/grid_y; // sub height
var c = document.createElement("canvas");
c.setAttribute("width",w);
c.setAttribute("Height",h);
var element = document.getElementById("insertcanvashere_freq");
element.appendChild(c);
var ctx = c.getContext("2d");
ctx.font="20px Georgia";
ctx.save();
ctx.setLineDash([1, 3]);/*dashes are 5px and spaces are 3px*/
ctx.strokeStyle="blue";
var skip_x = w/(frequency*2);
var y, max_y;
// draw x grid
var half = Math.floor(grid_y/2);
var skip_y = volt_size / half;
var voltage = 1;
for (var i=0; i<grid_y; i++) {
y=(-half+i)*(h/grid_y)+h/2;
if (i==0) {
max_y=y;
}
ctx.beginPath();
ctx.moveTo(0,y);
ctx.lineTo(w,y);
ctx.stroke();
// write voltages
if (i==half-1) {
//ctx.fillText(skip_y.toString() + 'V',skip_x * 2,y - 5);
}
}
min_y=y;
// draw y grid
var period=total_time / frequency;
var total_period=0;
var half = Math.ceil(grid_x/2);
for (var x=skip_x; x<w; x+=skip_x) {
ctx.beginPath();
ctx.moveTo(x,0);
ctx.lineTo(x,h);
ctx.stroke();
// write time period
total_period+=period;
half--;
if (half==0) {
round_period = Math.round(total_period * 1000) / 1000;
ctx.fillText(round_period.toString() + 'ms',x,min_y + 10);
// draw different color for the period
ctx.save();
ctx.setLineDash([10, 6]);
ctx.strokeStyle="black";
ctx.stroke();
ctx.restore();
}
}
ctx.restore();
// draw a sine wave
ctx.beginPath();
ctx.moveTo(0,h/2);
for (var x=1; x<w; x++) {
y=hs/2*Math.sin(2*Math.PI*frequency*x/w)+h/2;
ctx.lineTo(x,y);
}
ctx.strokeStyle="red";
ctx.stroke();
2.3 ניתוח גל סינוס - מתח אפקטיבי
var total_time = {total_time};
var frequency = {frequency};
var volt_size = {volt_size};
var w=500;
var h=360;
var grid_y=Math.round(frequency+4);
var grid_x=frequency*2-1;
var hs=h*(grid_y-1)/grid_y; // sub height
var c = document.createElement("canvas");
c.setAttribute("width",w);
c.setAttribute("Height",h);
var element = document.getElementById("insertcanvashere_volt");
element.appendChild(c);
var ctx = c.getContext("2d");
ctx.font="20px Georgia";
ctx.save();
ctx.setLineDash([1, 3]);/*dashes are 5px and spaces are 3px*/
ctx.strokeStyle="blue";
var skip_x = w/(frequency*2);
var y, max_y;
// draw x grid
var half = Math.floor(grid_y/2);
var skip_y = volt_size / ((grid_y-1)/2);
var voltage = 1;
for (var i=0; i<grid_y; i++) {
y=(-half+i)*(h/grid_y)+h/2;
if (i==0) {
max_y=y;
}
ctx.beginPath();
ctx.moveTo(0,y);
ctx.lineTo(w,y);
ctx.stroke();
// write voltages
if (i==half) {
ctx.fillText('0',w/2,y - 5);
// draw different color for the period
ctx.save();
ctx.setLineDash([1, 0]);
ctx.strokeStyle="black";
ctx.stroke();
ctx.restore();
}
if (i==half-1) {
round_volt = Math.round(skip_y * 1000) / 1000;
ctx.fillText(round_volt.toString() + 'V',w/2,y - 5);
// draw different color for the period
ctx.save();
ctx.setLineDash([10, 6]);
ctx.strokeStyle="black";
ctx.stroke();
ctx.restore();
}
}
min_y=y;
// draw y grid
var period=total_time / frequency;
var total_period=0;
var half = Math.ceil(grid_x/2);
for (var x=skip_x; x<w; x+=skip_x) {
ctx.beginPath();
ctx.moveTo(x,0);
ctx.lineTo(x,h);
ctx.stroke();
// write time period
total_period+=period;
half--;
if (half==1000) {
round_period = Math.round(total_period * 1000) / 1000;
ctx.fillText(round_period.toString() + 'ms',x,min_y + 10);
// draw different color for the period
ctx.save();
ctx.setLineDash([10, 6]);
ctx.strokeStyle="black";
ctx.stroke();
ctx.restore();
}
}
ctx.restore();
// draw a sine wave
ctx.beginPath();
ctx.moveTo(0,h/2);
for (var x=1; x<w; x++) {
y=hs/2*Math.sin(2*Math.PI*frequency*x/w)+h/2;
ctx.lineTo(x,y);
}
ctx.strokeStyle="red";
ctx.stroke();
3 מאקסימה
3.1 שליפת תדר בסריקה ליניארית
/* CAS Maxima to be used in moodle as a STACK question type */
/* Create a random frequency response of an RC circuit - either HPF or LPF */
rnd : rand(1.0);
/* high pass = 1, low pass = 0 */
high_pass : round(rnd);
is_high_pass: if high_pass = 1 then true else false;
f_cutoff: rand(41)+10; /* cutoff between 10 and 50 Hz */
/* time constant */
tau: 1/(2*pi*f_cutoff);
/* create omega RC (depends on f) */
wrc: 2*pi*f*tau;
/* Get the denominator */
denom: if high_pass=0 then 1 else wrc;
/* maximum voltage between 2 to 12 in steps of 0.5 */
volt_size : round((10*rnd+2)*2)/2;
/* absolute value of filter's output voltage */
v_out: volt_size*abs(denom/sqrt(1+wrc^2));
/* cutoff voltage */
v_cutoff: volt_size / sqrt(2);
f0: rand(10)+1; /* start frequency between 1 and 10 */
fT: rand(31)+70; /* stop frequency between 70 and 100 */
T: rand(6)+5; /* total time of frequency sweep between 5 and 10 seconds*/
tC: (f_cutoff-f0)*T/(fT-f0); /* cursor time position */
/* calculate f according to time */
f: (fT-f0)*t/T+f0;
/* do the plot */
p: plot(v_out, [t, 0, T], [ylabel, "Voltage [V]"], [xlabel, "Timing [S]"], grid2d, [y, 0, volt_size], [xtics, tC, T, T], [ytics, 0, v_cutoff, v_cutoff], [label, ["0", 0, 0], [string(T), T, 0]]);
3.2 תדר תהודה
/* CAS Maxima to be used in moodle as a STACK question type */
/* Create a random frequency response of a parallel RLC circuit */
Q: (rand(20)+6)/2; /* Q between 3 and 12.5 */
C: 0.01e-6; /* constant capacitor of 0.01uF */
w0: (rand(21)+10)*500 /* natural frequency between 5k and 15k */
L: 1/(C*w0^2) /* w0=1/sqrt(L*C) */
R: Q*sqrt(L/C); /* calculate R from Q */
alpha: 1/(2*R*C); /* damping ratio */
Vin: (rand(4)+1)/2; /* Input voltage between 0.5 and 2 */
Vout: Vin*R/sqrt((R-w^2*R*L*C)^2+(w*L)^2);
p: plot(Vout, [w, 0, w0+5*alpha], grid2d, [xlabel, "Angular frequency [rad/sec]"],[ylabel, "Voltage [V]"] );