I’ve been having trouble getting the DC Motor block’s results for current and rotational velocity to match what I obtain when using coupled differential equations in an ODE45 script. The results are very far apart and I’m not sure what is causing this.
The systems of equations I’ve implemented in ODE45 is standard in the DC motor literature, and is as follows:

Where is the current (), is the armature inductance (), is the armature resistance (), is the motor electrical constant (), is the voltage of the constant voltage source (), is the armature’s mass moment of inertia (), is the motor’s torque constant (), and is the motor damping ().

However, I’ve noticed that in the DC Motor block documentation, the equation for torque across the motor is given as (where is the damping constant). When I add the first term into my differential equation model, the results of ODE45 and Simscape Electrical become very close. I assume that this is the set of equations that Simscape is using, but I don’t like to use equations without understanding them. Where does this added term come from? I don’t believe I’ve seen it in the literature.

I’ve attached plots showing current and rotational velocity of the motor for the three models, so that you can see the discrepancies. (Note that there is zero damping in any of the models, yet only the classical differential equation model shows a rising angular velocity — the other models eventually reach some steady-state non-zero angular velocity.)

Any insight into what is happening here would be appreciated.I’ve been having trouble getting the DC Motor block’s results for current and rotational velocity to match what I obtain when using coupled differential equations in an ODE45 script. The results are very far apart and I’m not sure what is causing this.
The systems of equations I’ve implemented in ODE45 is standard in the DC motor literature, and is as follows:

Where is the current (), is the armature inductance (), is the armature resistance (), is the motor electrical constant (), is the voltage of the constant voltage source (), is the armature’s mass moment of inertia (), is the motor’s torque constant (), and is the motor damping ().

However, I’ve noticed that in the DC Motor block documentation, the equation for torque across the motor is given as (where is the damping constant). When I add the first term into my differential equation model, the results of ODE45 and Simscape Electrical become very close. I assume that this is the set of equations that Simscape is using, but I don’t like to use equations without understanding them. Where does this added term come from? I don’t believe I’ve seen it in the literature.

I’ve attached plots showing current and rotational velocity of the motor for the three models, so that you can see the discrepancies. (Note that there is zero damping in any of the models, yet only the classical differential equation model shows a rising angular velocity — the other models eventually reach some steady-state non-zero angular velocity.)

Any insight into what is happening here would be appreciated. I’ve been having trouble getting the DC Motor block’s results for current and rotational velocity to match what I obtain when using coupled differential equations in an ODE45 script. The results are very far apart and I’m not sure what is causing this.
The systems of equations I’ve implemented in ODE45 is standard in the DC motor literature, and is as follows:

Where is the current (), is the armature inductance (), is the armature resistance (), is the motor electrical constant (), is the voltage of the constant voltage source (), is the armature’s mass moment of inertia (), is the motor’s torque constant (), and is the motor damping ().

However, I’ve noticed that in the DC Motor block documentation, the equation for torque across the motor is given as (where is the damping constant). When I add the first term into my differential equation model, the results of ODE45 and Simscape Electrical become very close. I assume that this is the set of equations that Simscape is using, but I don’t like to use equations without understanding them. Where does this added term come from? I don’t believe I’ve seen it in the literature.

I’ve attached plots showing current and rotational velocity of the motor for the three models, so that you can see the discrepancies. (Note that there is zero damping in any of the models, yet only the classical differential equation model shows a rising angular velocity — the other models eventually reach some steady-state non-zero angular velocity.)

Any insight into what is happening here would be appreciated. simscape, differential equations MATLAB Answers — New Questions

​

I have error message while trying to link Matlab with my Aspen Plus simulation:
>> CHLC
Error using COM.Apwn_Document_37_0/invoke
Invoke Error, Dispatch Exception:
Source: Aspen Plus 37.0 OLE Services
Description: Unable to open file.

Error in CHLC (line 8)
Aspen.invoke(‘InitFromFile’,[mess.Name ” Simulation_Name ‘.apwz’]);
This is the code im using:
%% Linking
Aspen = actxserver(‘Apwn.Document.37.0’);
[stat,mess]=fileattrib; % get attributes of folder (Necessary to establish the location of the simulation)
Simulation_Name = ‘CL_SCO2’;% Aspeen Plus Simulation Name
Aspen.invoke(‘InitFromArchive2’,[mess.Name ” Simulation_Name ‘.bkp’]);
Aspen.Visible = 1; % 1 —> Aspen is Visible; 0 —> Aspen is open but not visible
Aspen.SuppressDialogs = 1; % Suppress windows dialogs.
Aspen.Engine.Run2(1); % Run the simulation
while Aspen.Engine.IsRunning == 1 % 1 –> If Aspen is running; 0 —> If Aspen stop.
pause(0.5);
end
Matlab version is R2022a (9.12.0.1884302) 64-bit and Aspen plus V11 (37.0.0395)
Already try saving Aspen simulation as .bkp and .apwz.
Yes, two files ( .m matlab code and .awpz / .bkp simulation) are in the same folder.I have error message while trying to link Matlab with my Aspen Plus simulation:
>> CHLC
Error using COM.Apwn_Document_37_0/invoke
Invoke Error, Dispatch Exception:
Source: Aspen Plus 37.0 OLE Services
Description: Unable to open file.

Error in CHLC (line 8)
Aspen.invoke(‘InitFromFile’,[mess.Name ” Simulation_Name ‘.apwz’]);
This is the code im using:
%% Linking
Aspen = actxserver(‘Apwn.Document.37.0’);
[stat,mess]=fileattrib; % get attributes of folder (Necessary to establish the location of the simulation)
Simulation_Name = ‘CL_SCO2’;% Aspeen Plus Simulation Name
Aspen.invoke(‘InitFromArchive2’,[mess.Name ” Simulation_Name ‘.bkp’]);
Aspen.Visible = 1; % 1 —> Aspen is Visible; 0 —> Aspen is open but not visible
Aspen.SuppressDialogs = 1; % Suppress windows dialogs.
Aspen.Engine.Run2(1); % Run the simulation
while Aspen.Engine.IsRunning == 1 % 1 –> If Aspen is running; 0 —> If Aspen stop.
pause(0.5);
end
Matlab version is R2022a (9.12.0.1884302) 64-bit and Aspen plus V11 (37.0.0395)
Already try saving Aspen simulation as .bkp and .apwz.
Yes, two files ( .m matlab code and .awpz / .bkp simulation) are in the same folder. I have error message while trying to link Matlab with my Aspen Plus simulation:
>> CHLC
Error using COM.Apwn_Document_37_0/invoke
Invoke Error, Dispatch Exception:
Source: Aspen Plus 37.0 OLE Services
Description: Unable to open file.

Error in CHLC (line 8)
Aspen.invoke(‘InitFromFile’,[mess.Name ” Simulation_Name ‘.apwz’]);
This is the code im using:
%% Linking
Aspen = actxserver(‘Apwn.Document.37.0’);
[stat,mess]=fileattrib; % get attributes of folder (Necessary to establish the location of the simulation)
Simulation_Name = ‘CL_SCO2’;% Aspeen Plus Simulation Name
Aspen.invoke(‘InitFromArchive2’,[mess.Name ” Simulation_Name ‘.bkp’]);
Aspen.Visible = 1; % 1 —> Aspen is Visible; 0 —> Aspen is open but not visible
Aspen.SuppressDialogs = 1; % Suppress windows dialogs.
Aspen.Engine.Run2(1); % Run the simulation
while Aspen.Engine.IsRunning == 1 % 1 –> If Aspen is running; 0 —> If Aspen stop.
pause(0.5);
end
Matlab version is R2022a (9.12.0.1884302) 64-bit and Aspen plus V11 (37.0.0395)
Already try saving Aspen simulation as .bkp and .apwz.
Yes, two files ( .m matlab code and .awpz / .bkp simulation) are in the same folder. aspen puls, matlab, link, error MATLAB Answers — New Questions

​

how do i model power increment (+udeg, -udeg) limits using matlab function?(sign function)
i am working on power system having constraints.
Thanks in advancehow do i model power increment (+udeg, -udeg) limits using matlab function?(sign function)
i am working on power system having constraints.
Thanks in advance how do i model power increment (+udeg, -udeg) limits using matlab function?(sign function)
i am working on power system having constraints.
Thanks in advance sign function MATLAB Answers — New Questions

​

Hello,
I am working with the IEEE 34-bus distribution network in MATLAB-Simulink. When running the load flow analysis after setting up the buses, I encountered the following error message:
"Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.218374e-17."
Given that this is an unbalanced system, could you please advise on the potential causes of this issue?How to solve it? Additionally, are there alternative methods or tools within MATLAB-Simulink to handle load flow analysis for such systems?
Thank you for your assistance.
Best regards,Hello,
I am working with the IEEE 34-bus distribution network in MATLAB-Simulink. When running the load flow analysis after setting up the buses, I encountered the following error message:
"Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.218374e-17."
Given that this is an unbalanced system, could you please advise on the potential causes of this issue?How to solve it? Additionally, are there alternative methods or tools within MATLAB-Simulink to handle load flow analysis for such systems?
Thank you for your assistance.
Best regards, Hello,
I am working with the IEEE 34-bus distribution network in MATLAB-Simulink. When running the load flow analysis after setting up the buses, I encountered the following error message:
"Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.218374e-17."
Given that this is an unbalanced system, could you please advise on the potential causes of this issue?How to solve it? Additionally, are there alternative methods or tools within MATLAB-Simulink to handle load flow analysis for such systems?
Thank you for your assistance.
Best regards, ieee 34, load flow, matlab, simulink, unbalanced s MATLAB Answers — New Questions

​

I would like to speed up the following code which solves a discrete dynamic programming problem using the method of successive approximations, as described e.g. in Bertsekas.
The algorithm is made of two steps. In step 1, I precompute the payoff array R(a’,a,z) where a’ is the action and (a,z) are the states. In step 2 I compute the value function using the method of successive approximations: I guess V0, then I compute an updated V1 and finally I check if ||V1-V0|| is less than a tolerance level. If it is, I stop, otherwise I set V0=V1 and go on.
I profiled the code (see below a MWE) and the two most time-consuming lines are the following ones:
(1) RHS = Ret+beta*permute(EV,[1,3,2]);
(2) [max_val,max_ind] = max(RHS,[],1);
Line (1) takes up 63% of the total running time, line (2) takes 32%. As you can see I have already vectorized all loops.
I would be very grateful for any suggestion. I post below a MWE.
%% Solve income fluctuation problem CPU
clear;clc;close all
%% Economic parameters
sigma = 2;
r = 0.03;
beta = 0.96;
PZ = [0.60 0.40;
0.05 0.95];
z_grid = [0.5 1.0]’;
n_z = length(z_grid);
b = 0; %lower bound for asset holdings a
grid_max = 4;
n_a = 500; % IN PRACTICE THIS IS EQUAL TO 5000-10000
R = 1+r;
a_grid = linspace(-b,grid_max,n_a)’;
if sigma==1
fun_u = @(c) log(c);
else
fun_u = @(c) c.^(1-sigma)/(1-sigma);
end
%% Computational parameters
verbose = 0;
tiny = 1e-8; %very small positive number
tol = 1e-6; %tolerance for VFI and TI
max_iter = 500; %maximum num. of iterations for both VFI and TI
%% Start timing
tic
%% STEP 1- Precompute current payoff array R(a’,a,z)
a_tomorrow = a_grid; %(a’,1,1)
a_today = a_grid’; %(1,a,1)
z_today = shiftdim(z_grid,-2); %(1,1,z)
cons = (1+r)*a_today+z_today-a_tomorrow;
Ret = fun_u(cons); %size: [n_a,n_a,n_z]
Ret(cons<=0) = -inf;
%% STEP 2 – Value function iteration
iter = 1;
err = tol+1;
V0 = zeros(n_a,n_z);
while err>tol && iter<=max_iter
EV = V0*PZ’; %(a’,z)
RHS = Ret+beta*permute(EV,[1,3,2]);
[max_val,max_ind] = max(RHS,[],1);
V1 = squeeze(max_val);
pol_ind_ap = squeeze(max_ind);
err = max(abs(V0(:)-V1(:)));
if verbose==1
fprintf(‘iter = %d, err = %f n’,iter,err)
end
iter = iter+1;
V0 = V1;
end
if err>tol
error(‘VFI did not converge!’)
else
fprintf(‘VFI converged after = %d iterations n’,iter)
end
pol_ap = a_grid(pol_ind_ap);
pol_c = (1+r)*a_grid+z_grid’-pol_ap;
%% End timing
toc
%% Figures
figure
plot(a_grid,pol_c(:,1),’linewidth’,2)
hold on
plot(a_grid,pol_c(:,2),’linewidth’,2)
legend(‘Low shock’,’High shock’,’Location’,’NorthWest’)
xlabel(‘asset level’)
ylabel(‘consumption’)
title(‘Consumption Policy Function’)
figure
plot(a_grid,a_grid,’–‘,’linewidth’,2)
hold on
plot(a_grid,pol_ap(:,1),’linewidth’,2)
hold on
plot(a_grid,pol_ap(:,2),’linewidth’,2)
legend(’45 line’,’Low shock’,’High shock’,’Location’,’NorthWest’)
xlabel(‘Current period assets’)
ylabel(‘Next-period assets’)
title(‘Assets Policy Function’)I would like to speed up the following code which solves a discrete dynamic programming problem using the method of successive approximations, as described e.g. in Bertsekas.
The algorithm is made of two steps. In step 1, I precompute the payoff array R(a’,a,z) where a’ is the action and (a,z) are the states. In step 2 I compute the value function using the method of successive approximations: I guess V0, then I compute an updated V1 and finally I check if ||V1-V0|| is less than a tolerance level. If it is, I stop, otherwise I set V0=V1 and go on.
I profiled the code (see below a MWE) and the two most time-consuming lines are the following ones:
(1) RHS = Ret+beta*permute(EV,[1,3,2]);
(2) [max_val,max_ind] = max(RHS,[],1);
Line (1) takes up 63% of the total running time, line (2) takes 32%. As you can see I have already vectorized all loops.
I would be very grateful for any suggestion. I post below a MWE.
%% Solve income fluctuation problem CPU
clear;clc;close all
%% Economic parameters
sigma = 2;
r = 0.03;
beta = 0.96;
PZ = [0.60 0.40;
0.05 0.95];
z_grid = [0.5 1.0]’;
n_z = length(z_grid);
b = 0; %lower bound for asset holdings a
grid_max = 4;
n_a = 500; % IN PRACTICE THIS IS EQUAL TO 5000-10000
R = 1+r;
a_grid = linspace(-b,grid_max,n_a)’;
if sigma==1
fun_u = @(c) log(c);
else
fun_u = @(c) c.^(1-sigma)/(1-sigma);
end
%% Computational parameters
verbose = 0;
tiny = 1e-8; %very small positive number
tol = 1e-6; %tolerance for VFI and TI
max_iter = 500; %maximum num. of iterations for both VFI and TI
%% Start timing
tic
%% STEP 1- Precompute current payoff array R(a’,a,z)
a_tomorrow = a_grid; %(a’,1,1)
a_today = a_grid’; %(1,a,1)
z_today = shiftdim(z_grid,-2); %(1,1,z)
cons = (1+r)*a_today+z_today-a_tomorrow;
Ret = fun_u(cons); %size: [n_a,n_a,n_z]
Ret(cons<=0) = -inf;
%% STEP 2 – Value function iteration
iter = 1;
err = tol+1;
V0 = zeros(n_a,n_z);
while err>tol && iter<=max_iter
EV = V0*PZ’; %(a’,z)
RHS = Ret+beta*permute(EV,[1,3,2]);
[max_val,max_ind] = max(RHS,[],1);
V1 = squeeze(max_val);
pol_ind_ap = squeeze(max_ind);
err = max(abs(V0(:)-V1(:)));
if verbose==1
fprintf(‘iter = %d, err = %f n’,iter,err)
end
iter = iter+1;
V0 = V1;
end
if err>tol
error(‘VFI did not converge!’)
else
fprintf(‘VFI converged after = %d iterations n’,iter)
end
pol_ap = a_grid(pol_ind_ap);
pol_c = (1+r)*a_grid+z_grid’-pol_ap;
%% End timing
toc
%% Figures
figure
plot(a_grid,pol_c(:,1),’linewidth’,2)
hold on
plot(a_grid,pol_c(:,2),’linewidth’,2)
legend(‘Low shock’,’High shock’,’Location’,’NorthWest’)
xlabel(‘asset level’)
ylabel(‘consumption’)
title(‘Consumption Policy Function’)
figure
plot(a_grid,a_grid,’–‘,’linewidth’,2)
hold on
plot(a_grid,pol_ap(:,1),’linewidth’,2)
hold on
plot(a_grid,pol_ap(:,2),’linewidth’,2)
legend(’45 line’,’Low shock’,’High shock’,’Location’,’NorthWest’)
xlabel(‘Current period assets’)
ylabel(‘Next-period assets’)
title(‘Assets Policy Function’) I would like to speed up the following code which solves a discrete dynamic programming problem using the method of successive approximations, as described e.g. in Bertsekas.
The algorithm is made of two steps. In step 1, I precompute the payoff array R(a’,a,z) where a’ is the action and (a,z) are the states. In step 2 I compute the value function using the method of successive approximations: I guess V0, then I compute an updated V1 and finally I check if ||V1-V0|| is less than a tolerance level. If it is, I stop, otherwise I set V0=V1 and go on.
I profiled the code (see below a MWE) and the two most time-consuming lines are the following ones:
(1) RHS = Ret+beta*permute(EV,[1,3,2]);
(2) [max_val,max_ind] = max(RHS,[],1);
Line (1) takes up 63% of the total running time, line (2) takes 32%. As you can see I have already vectorized all loops.
I would be very grateful for any suggestion. I post below a MWE.
%% Solve income fluctuation problem CPU
clear;clc;close all
%% Economic parameters
sigma = 2;
r = 0.03;
beta = 0.96;
PZ = [0.60 0.40;
0.05 0.95];
z_grid = [0.5 1.0]’;
n_z = length(z_grid);
b = 0; %lower bound for asset holdings a
grid_max = 4;
n_a = 500; % IN PRACTICE THIS IS EQUAL TO 5000-10000
R = 1+r;
a_grid = linspace(-b,grid_max,n_a)’;
if sigma==1
fun_u = @(c) log(c);
else
fun_u = @(c) c.^(1-sigma)/(1-sigma);
end
%% Computational parameters
verbose = 0;
tiny = 1e-8; %very small positive number
tol = 1e-6; %tolerance for VFI and TI
max_iter = 500; %maximum num. of iterations for both VFI and TI
%% Start timing
tic
%% STEP 1- Precompute current payoff array R(a’,a,z)
a_tomorrow = a_grid; %(a’,1,1)
a_today = a_grid’; %(1,a,1)
z_today = shiftdim(z_grid,-2); %(1,1,z)
cons = (1+r)*a_today+z_today-a_tomorrow;
Ret = fun_u(cons); %size: [n_a,n_a,n_z]
Ret(cons<=0) = -inf;
%% STEP 2 – Value function iteration
iter = 1;
err = tol+1;
V0 = zeros(n_a,n_z);
while err>tol && iter<=max_iter
EV = V0*PZ’; %(a’,z)
RHS = Ret+beta*permute(EV,[1,3,2]);
[max_val,max_ind] = max(RHS,[],1);
V1 = squeeze(max_val);
pol_ind_ap = squeeze(max_ind);
err = max(abs(V0(:)-V1(:)));
if verbose==1
fprintf(‘iter = %d, err = %f n’,iter,err)
end
iter = iter+1;
V0 = V1;
end
if err>tol
error(‘VFI did not converge!’)
else
fprintf(‘VFI converged after = %d iterations n’,iter)
end
pol_ap = a_grid(pol_ind_ap);
pol_c = (1+r)*a_grid+z_grid’-pol_ap;
%% End timing
toc
%% Figures
figure
plot(a_grid,pol_c(:,1),’linewidth’,2)
hold on
plot(a_grid,pol_c(:,2),’linewidth’,2)
legend(‘Low shock’,’High shock’,’Location’,’NorthWest’)
xlabel(‘asset level’)
ylabel(‘consumption’)
title(‘Consumption Policy Function’)
figure
plot(a_grid,a_grid,’–‘,’linewidth’,2)
hold on
plot(a_grid,pol_ap(:,1),’linewidth’,2)
hold on
plot(a_grid,pol_ap(:,2),’linewidth’,2)
legend(’45 line’,’Low shock’,’High shock’,’Location’,’NorthWest’)
xlabel(‘Current period assets’)
ylabel(‘Next-period assets’)
title(‘Assets Policy Function’) vectorized, performance MATLAB Answers — New Questions

​

Hello, I am working on this matlab script for my final paper and simulating the transmission of 1 pulse. First it gets shaped as a Raised Cossine pulse with a filter, then the transmitted Tx signal goes through a low-pass butterworth filter that acts like a channel introducing distortion to the pulse, then the received Rx signal comes to a Zero Forcing Equalizer filter that should correct the low-pass filter distortion.
The problem I’m facing is when I try to plot the pulse, here’s the code:
clear;clf;clc;close all;

% This script performs a transmission of 1 pulse shaped with Raised Cossine
% that receives a distortion from a butterworth low-pass filter then
% the pulse gets an equalization from a Zero Forcing filter

% Building RC pulse shaping filter
rf = 0.5; % Roll-off factor
span = 16; % Pulse truncation factor
sps = 10; % Samples per pulse
p_rc = rcosdesign(rf,span,sps,"normal"); % Designing filter
p_rc_norm = p_rc/max(p_rc); % Normalizing
h_span = span/2; % Half symbol span

% Setting Data parameters
L = 1; % Total Data pulses
Rb = 1000; % Data rate in bps
% For Rb = 1000, that means that symbol time -> Tb = 1 milisecond
Fs = Rb * sps; % Sampling frequency
x = 1; % Data signal
% Note: Not starting at zero because of filter group delay
tx = 1000*(h_span:h_span)/Rb; % x Time vector sampled
% at symbol rate in milliseconds -> Data

% Applying RC filter
x_rc = upfirdn(x, p_rc_norm, sps); % Upsampling with upfirdn
% Note: By applying filtering, remember that the peak response of the
% filter gets delayed by the group delay of the filter, so we need to
% get rid of transients so that the pulses can be overlapping when plotting
trans_1 = sps*h_span; % Number of points for transient before peak operation
trans_2 = sps*(h_span-1); % Number of points for transient after peak operation
tx_rc = 1000*(0:L*sps+trans_1+trans_2)/Fs; % x_rc Time vector sampled at
% sampling frequency in milliseconds -> Tx signal

% Building and applying channel distortion via butterworth low-pass filter
Wn_butter = 0.1; % Normalized cutoff frequency
n_butter = 4; % Butterworth nth-order
[b_chan,a_chan] = butter(n_butter, Wn_butter); % Low-pass filter coefficients
x_rc_ch = filtfilt(b_chan,a_chan,x_rc); % Used filtfilt() so that it doesn’t
% get any group delay when filtering

% Setting Zero Forcing Equalizer filter parameters
N2_plus_1 = 9; % Choose ZF EQ Filter Order = 2N+1 % Note: Odd number only
N2 = N2_plus_1 – 1;
% If ZF EQ Filter Order = 5, that means that N = 2, that also
% means that the matrix would be 5×5, because the row and column
% should be from Pr[0] to Pr[2N] = Pr[4] and from Pr[0] to Pr[-2N]= Pr[-4]

% Building and applying ZFE filter
zf_kTb_peak = max(x_rc_ch); % Identify pulse peak
zf_kTb_peak_pos = find(x_rc_ch==zf_kTb_peak); % Identify pulse peak position
Pr_minus_N2 = zf_kTb_peak_pos-(sps*N2); % Range points to feed the ZF equalizer: Pr[0] to Pr[-2N]
Pr_plus_N2 = zf_kTb_peak_pos+(sps*N2); % Range points to feed the ZF equalizer: Pr[0] to Pr[2N]
zf_kTb = x_rc_ch(Pr_minus_N2:sps:Pr_plus_N2); % Scrolling Rx signal
% to get kTb points (integer multiples"k" of symbol time"Tb")
zf_kTb = zf_kTb’;
zf_Pr_row = zf_kTb(1:round(length(zf_kTb)/2));
zf_Pr_row = zf_Pr_row(end:-1:1); % Getting Pr[0] to Pr[-2N] vector
zf_Pr_col = zf_kTb(round(length(zf_kTb)/2):end); % Getting Pr[0] to Pr[2N] vector
Pr = toeplitz(zf_Pr_col, zf_Pr_row); % With kTb points we can build Toeplitz Matrix
Po = zeros(round(length(zf_kTb)/2), 1); % Creating Po vector to force neighbour pulses to 0
Po(round(length(Po)/2)) = 1;
Ck_zf = inv(Pr)*Po; % Calculating ZFE filter coefficients

% Applying ZFE filter to Rx kTb points:
zf_kTb_ZFE = conv(zf_kTb, Ck_zf, ‘same’);
% Plot 1
figure(1)
stem(zf_kTb, ‘color’, [0.8500 0.3250 0.0980]);
hold on
stem(zf_kTb_ZFE, ‘color’, [0.9290 0.6940 0.1250]);
hold off;
xlabel(‘Symbol spaced points (kTb)’); ylabel(‘Magnitude’);
legend(‘Rx points’,’ZFE points’,’Location’,’southeast’);

% Applying ZFE filter to Rx signal
x_rc_ch_zf = conv(x_rc_ch, Ck_zf, ‘same’);
% Plot 2
figure(2)
stem(tx, x, ‘k’);
hold on;
plot(tx_rc, x_rc, ‘color’, [0 0.4470 0.7410]);
plot(tx_rc, x_rc_ch, ‘color’, [0.8500 0.3250 0.0980]);
plot(tx_rc, x_rc_ch_zf, ‘color’, [0.9290 0.6940 0.1250]);
hold off;
xlabel(‘Time (ms)’); ylabel(‘Magnitude’);
legend(‘Data’,’Tx Signal’,’Rx Signal’,’ZFE Signal’,’Location’,’southeast’);
For the first figure, the equalyzer is working like it should, it is forcing the symbol spaced points(kTb) to zero.
But for the second figure, when the whole discretized pulse from the Rx signal goes through the same ZFE filter, it doesn’t get equalyzed back to something like the Tx signal, it just gets a little bit bigger on the peak.
So my doubt is, can we equalyze the Rx pulse with this type of equalization and plot the equalized pulse to something close to Tx signal ?Or does it just works for the points of symbol time ?Hello, I am working on this matlab script for my final paper and simulating the transmission of 1 pulse. First it gets shaped as a Raised Cossine pulse with a filter, then the transmitted Tx signal goes through a low-pass butterworth filter that acts like a channel introducing distortion to the pulse, then the received Rx signal comes to a Zero Forcing Equalizer filter that should correct the low-pass filter distortion.
The problem I’m facing is when I try to plot the pulse, here’s the code:
clear;clf;clc;close all;

% This script performs a transmission of 1 pulse shaped with Raised Cossine
% that receives a distortion from a butterworth low-pass filter then
% the pulse gets an equalization from a Zero Forcing filter

% Building RC pulse shaping filter
rf = 0.5; % Roll-off factor
span = 16; % Pulse truncation factor
sps = 10; % Samples per pulse
p_rc = rcosdesign(rf,span,sps,"normal"); % Designing filter
p_rc_norm = p_rc/max(p_rc); % Normalizing
h_span = span/2; % Half symbol span

% Setting Data parameters
L = 1; % Total Data pulses
Rb = 1000; % Data rate in bps
% For Rb = 1000, that means that symbol time -> Tb = 1 milisecond
Fs = Rb * sps; % Sampling frequency
x = 1; % Data signal
% Note: Not starting at zero because of filter group delay
tx = 1000*(h_span:h_span)/Rb; % x Time vector sampled
% at symbol rate in milliseconds -> Data

% Applying RC filter
x_rc = upfirdn(x, p_rc_norm, sps); % Upsampling with upfirdn
% Note: By applying filtering, remember that the peak response of the
% filter gets delayed by the group delay of the filter, so we need to
% get rid of transients so that the pulses can be overlapping when plotting
trans_1 = sps*h_span; % Number of points for transient before peak operation
trans_2 = sps*(h_span-1); % Number of points for transient after peak operation
tx_rc = 1000*(0:L*sps+trans_1+trans_2)/Fs; % x_rc Time vector sampled at
% sampling frequency in milliseconds -> Tx signal

% Building and applying channel distortion via butterworth low-pass filter
Wn_butter = 0.1; % Normalized cutoff frequency
n_butter = 4; % Butterworth nth-order
[b_chan,a_chan] = butter(n_butter, Wn_butter); % Low-pass filter coefficients
x_rc_ch = filtfilt(b_chan,a_chan,x_rc); % Used filtfilt() so that it doesn’t
% get any group delay when filtering

% Setting Zero Forcing Equalizer filter parameters
N2_plus_1 = 9; % Choose ZF EQ Filter Order = 2N+1 % Note: Odd number only
N2 = N2_plus_1 – 1;
% If ZF EQ Filter Order = 5, that means that N = 2, that also
% means that the matrix would be 5×5, because the row and column
% should be from Pr[0] to Pr[2N] = Pr[4] and from Pr[0] to Pr[-2N]= Pr[-4]

% Building and applying ZFE filter
zf_kTb_peak = max(x_rc_ch); % Identify pulse peak
zf_kTb_peak_pos = find(x_rc_ch==zf_kTb_peak); % Identify pulse peak position
Pr_minus_N2 = zf_kTb_peak_pos-(sps*N2); % Range points to feed the ZF equalizer: Pr[0] to Pr[-2N]
Pr_plus_N2 = zf_kTb_peak_pos+(sps*N2); % Range points to feed the ZF equalizer: Pr[0] to Pr[2N]
zf_kTb = x_rc_ch(Pr_minus_N2:sps:Pr_plus_N2); % Scrolling Rx signal
% to get kTb points (integer multiples"k" of symbol time"Tb")
zf_kTb = zf_kTb’;
zf_Pr_row = zf_kTb(1:round(length(zf_kTb)/2));
zf_Pr_row = zf_Pr_row(end:-1:1); % Getting Pr[0] to Pr[-2N] vector
zf_Pr_col = zf_kTb(round(length(zf_kTb)/2):end); % Getting Pr[0] to Pr[2N] vector
Pr = toeplitz(zf_Pr_col, zf_Pr_row); % With kTb points we can build Toeplitz Matrix
Po = zeros(round(length(zf_kTb)/2), 1); % Creating Po vector to force neighbour pulses to 0
Po(round(length(Po)/2)) = 1;
Ck_zf = inv(Pr)*Po; % Calculating ZFE filter coefficients

% Applying ZFE filter to Rx kTb points:
zf_kTb_ZFE = conv(zf_kTb, Ck_zf, ‘same’);
% Plot 1
figure(1)
stem(zf_kTb, ‘color’, [0.8500 0.3250 0.0980]);
hold on
stem(zf_kTb_ZFE, ‘color’, [0.9290 0.6940 0.1250]);
hold off;
xlabel(‘Symbol spaced points (kTb)’); ylabel(‘Magnitude’);
legend(‘Rx points’,’ZFE points’,’Location’,’southeast’);

% Applying ZFE filter to Rx signal
x_rc_ch_zf = conv(x_rc_ch, Ck_zf, ‘same’);
% Plot 2
figure(2)
stem(tx, x, ‘k’);
hold on;
plot(tx_rc, x_rc, ‘color’, [0 0.4470 0.7410]);
plot(tx_rc, x_rc_ch, ‘color’, [0.8500 0.3250 0.0980]);
plot(tx_rc, x_rc_ch_zf, ‘color’, [0.9290 0.6940 0.1250]);
hold off;
xlabel(‘Time (ms)’); ylabel(‘Magnitude’);
legend(‘Data’,’Tx Signal’,’Rx Signal’,’ZFE Signal’,’Location’,’southeast’);
For the first figure, the equalyzer is working like it should, it is forcing the symbol spaced points(kTb) to zero.
But for the second figure, when the whole discretized pulse from the Rx signal goes through the same ZFE filter, it doesn’t get equalyzed back to something like the Tx signal, it just gets a little bit bigger on the peak.
So my doubt is, can we equalyze the Rx pulse with this type of equalization and plot the equalized pulse to something close to Tx signal ?Or does it just works for the points of symbol time ? Hello, I am working on this matlab script for my final paper and simulating the transmission of 1 pulse. First it gets shaped as a Raised Cossine pulse with a filter, then the transmitted Tx signal goes through a low-pass butterworth filter that acts like a channel introducing distortion to the pulse, then the received Rx signal comes to a Zero Forcing Equalizer filter that should correct the low-pass filter distortion.
The problem I’m facing is when I try to plot the pulse, here’s the code:
clear;clf;clc;close all;

% This script performs a transmission of 1 pulse shaped with Raised Cossine
% that receives a distortion from a butterworth low-pass filter then
% the pulse gets an equalization from a Zero Forcing filter

% Building RC pulse shaping filter
rf = 0.5; % Roll-off factor
span = 16; % Pulse truncation factor
sps = 10; % Samples per pulse
p_rc = rcosdesign(rf,span,sps,"normal"); % Designing filter
p_rc_norm = p_rc/max(p_rc); % Normalizing
h_span = span/2; % Half symbol span

% Setting Data parameters
L = 1; % Total Data pulses
Rb = 1000; % Data rate in bps
% For Rb = 1000, that means that symbol time -> Tb = 1 milisecond
Fs = Rb * sps; % Sampling frequency
x = 1; % Data signal
% Note: Not starting at zero because of filter group delay
tx = 1000*(h_span:h_span)/Rb; % x Time vector sampled
% at symbol rate in milliseconds -> Data

% Applying RC filter
x_rc = upfirdn(x, p_rc_norm, sps); % Upsampling with upfirdn
% Note: By applying filtering, remember that the peak response of the
% filter gets delayed by the group delay of the filter, so we need to
% get rid of transients so that the pulses can be overlapping when plotting
trans_1 = sps*h_span; % Number of points for transient before peak operation
trans_2 = sps*(h_span-1); % Number of points for transient after peak operation
tx_rc = 1000*(0:L*sps+trans_1+trans_2)/Fs; % x_rc Time vector sampled at
% sampling frequency in milliseconds -> Tx signal

% Building and applying channel distortion via butterworth low-pass filter
Wn_butter = 0.1; % Normalized cutoff frequency
n_butter = 4; % Butterworth nth-order
[b_chan,a_chan] = butter(n_butter, Wn_butter); % Low-pass filter coefficients
x_rc_ch = filtfilt(b_chan,a_chan,x_rc); % Used filtfilt() so that it doesn’t
% get any group delay when filtering

% Setting Zero Forcing Equalizer filter parameters
N2_plus_1 = 9; % Choose ZF EQ Filter Order = 2N+1 % Note: Odd number only
N2 = N2_plus_1 – 1;
% If ZF EQ Filter Order = 5, that means that N = 2, that also
% means that the matrix would be 5×5, because the row and column
% should be from Pr[0] to Pr[2N] = Pr[4] and from Pr[0] to Pr[-2N]= Pr[-4]

% Building and applying ZFE filter
zf_kTb_peak = max(x_rc_ch); % Identify pulse peak
zf_kTb_peak_pos = find(x_rc_ch==zf_kTb_peak); % Identify pulse peak position
Pr_minus_N2 = zf_kTb_peak_pos-(sps*N2); % Range points to feed the ZF equalizer: Pr[0] to Pr[-2N]
Pr_plus_N2 = zf_kTb_peak_pos+(sps*N2); % Range points to feed the ZF equalizer: Pr[0] to Pr[2N]
zf_kTb = x_rc_ch(Pr_minus_N2:sps:Pr_plus_N2); % Scrolling Rx signal
% to get kTb points (integer multiples"k" of symbol time"Tb")
zf_kTb = zf_kTb’;
zf_Pr_row = zf_kTb(1:round(length(zf_kTb)/2));
zf_Pr_row = zf_Pr_row(end:-1:1); % Getting Pr[0] to Pr[-2N] vector
zf_Pr_col = zf_kTb(round(length(zf_kTb)/2):end); % Getting Pr[0] to Pr[2N] vector
Pr = toeplitz(zf_Pr_col, zf_Pr_row); % With kTb points we can build Toeplitz Matrix
Po = zeros(round(length(zf_kTb)/2), 1); % Creating Po vector to force neighbour pulses to 0
Po(round(length(Po)/2)) = 1;
Ck_zf = inv(Pr)*Po; % Calculating ZFE filter coefficients

% Applying ZFE filter to Rx kTb points:
zf_kTb_ZFE = conv(zf_kTb, Ck_zf, ‘same’);
% Plot 1
figure(1)
stem(zf_kTb, ‘color’, [0.8500 0.3250 0.0980]);
hold on
stem(zf_kTb_ZFE, ‘color’, [0.9290 0.6940 0.1250]);
hold off;
xlabel(‘Symbol spaced points (kTb)’); ylabel(‘Magnitude’);
legend(‘Rx points’,’ZFE points’,’Location’,’southeast’);

% Applying ZFE filter to Rx signal
x_rc_ch_zf = conv(x_rc_ch, Ck_zf, ‘same’);
% Plot 2
figure(2)
stem(tx, x, ‘k’);
hold on;
plot(tx_rc, x_rc, ‘color’, [0 0.4470 0.7410]);
plot(tx_rc, x_rc_ch, ‘color’, [0.8500 0.3250 0.0980]);
plot(tx_rc, x_rc_ch_zf, ‘color’, [0.9290 0.6940 0.1250]);
hold off;
xlabel(‘Time (ms)’); ylabel(‘Magnitude’);
legend(‘Data’,’Tx Signal’,’Rx Signal’,’ZFE Signal’,’Location’,’southeast’);
For the first figure, the equalyzer is working like it should, it is forcing the symbol spaced points(kTb) to zero.
But for the second figure, when the whole discretized pulse from the Rx signal goes through the same ZFE filter, it doesn’t get equalyzed back to something like the Tx signal, it just gets a little bit bigger on the peak.
So my doubt is, can we equalyze the Rx pulse with this type of equalization and plot the equalized pulse to something close to Tx signal ?Or does it just works for the points of symbol time ? filter, plot, zero-forcing MATLAB Answers — New Questions

​