Using Help using filter function in Matlab

2 visualizaciones (últimos 30 días)
Telema Harry
Telema Harry el 27 de Nov. de 2020
Editada: Telema Harry el 29 de Nov. de 2020
clear all
G = @(u,t)(25-(5-(u)).^2);
u = 0;
%g0 = G(0.1,0);
phase = 0;
% Extremum Seeking Control Parameters
freq = 10*2*pi; % sample frequency
dt = 1/freq;
T = 10; % total period of simulation (in seconds)
A = .2; % amplitude
omega = 10*2*pi; % 10 Hz
phase = 0;
K = 5; % integration gain
%*****************************************************************
% Design of High Pass filter. the cut of frequency is 0.2
%a and b at the output variable from the butter function
b_order = 1; % The butterworth filter order
butter_freq = 1;
[b,a] = butter(b_order,butter_freq * 2 * dt, 'high');
%freqz(b,a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Intial condition
%ys = zeros(1,b_order+1);
uhat(1) = u(1);
for i = 1 : T/dt
t = (i-1)*dt;
yvals(i) = G(u,t);% This variable stores the values of G at every interation
ys(i) = yvals(i);
HPFnew = filter(b,a,ys(i)); % This is the value of the output signal after passing through the filter.
HPFnew_vals(i) = HPFnew;
xi = HPFnew.*sin(omega.*t + phase);
uhat = uhat + xi.*K.*dt;
u = uhat + A*sin(omega*t + phase);
uhats(i) = uhat;
uvals(i) = u;
end
t = dt:dt:T;
figure
plot(t, HPFnew_vals)
%figure
%fplot(@(u) (25-(5-(u)).^2),[0 T]) % Ploting the original function.
% The next step is to multiply the output cost function G with the filter
% I will call that variable rho
t = dt:dt:T;
% Plot of U against Time
figure
subplot(2,1,1)
%ax1 = nexttile;
%plot([1:length(uhat)],uhat,'r','lineWidth',1)
plot(t,uhats,'r','lineWidth',1)
hold on
%plot([1:i],uvals,'k','lineWidth',1)
plot(t,uvals,'k','lineWidth',1)
hold off
xlabel ('Time')
ylabel('U')
legend('Uhat','U')
grid on
subplot(2,1,2)
%ax2 = nexttile;
%plot([1:length(yvals)],yvals,'r','lineWidth',1)
plot(t,yvals,'b','lineWidth',1)
xlabel ('Time')
ylabel('G')
grid on
I am trying to pass my output signal yvals(i) = G(u,t) through a butter highpass filter. The code runs but my results is totally different from what I was expecting.
This is the result obtained.
Result Obtained
This is the result expected. I don't know what i'm doing wrong?
  3 comentarios
Telema Harry
Telema Harry el 27 de Nov. de 2020
I saw the code in a text book (Data-driven Science and Engineering by Steven Brunton). I tried developing mine because I do not understand some part of the code.
clear all
% Extremum Seeking Control Parameters
J = @(u,t)(25-(5-(u)).^2);
y0 = J(0,0); %
u = 0;
% Extremum Seeking Control Parameters
freq = 10*2*pi; % sample frequency
dt = 1/freq;
T = 10; % total period of simulation (in seconds)
A = .2; % amplitude
omega = 10*2*pi; % 10 Hz
phase = 0;
K = 5; % integration gain
% High pass filter (Butterworth filter)
butterorder=1;
butterfreq = 1; % in Hz for ’high’
[b,a] = butter(butterorder,butterfreq * 2 * dt,'high')
ys = zeros(1,butterorder+1)+y0;
HPF = zeros(1,butterorder+1);
uhat=u;
for i=1:T/dt
t = (i-1)*dt;
yvals(i)=J(u,t);
for k=1:butterorder
ys(k) = ys(k+1);
HPF(k) = HPF(k+1);
end
ys(butterorder+1) = yvals(i);
HPFnew = 0;
for k=1:butterorder+1
HPFnew = HPFnew + b(k)*ys(butterorder+2-k);
end
for k=2:butterorder+1
HPFnew = HPFnew - a(k)*HPF(butterorder+2-k);
end
HPF(butterorder+1) = HPFnew;
xi = HPFnew*sin(omega*t + phase);
uhat = uhat + xi*K*dt;
u = uhat + A*sin(omega*t + phase);
uhats(i) = uhat;
uvals(i) = u;
end
t = dt:dt:T;
% Plot of U against Time
subplot(2,1,1)
%ax1 = nexttile;
%plot([1:length(uhat)],uhat,'r','lineWidth',1)
plot(t,uhats,'r','lineWidth',1)
hold on
%plot([1:i],uvals,'k','lineWidth',1)
plot(t,uvals,'k','lineWidth',1)
hold off
xlabel ('Time')
ylabel('U')
legend('Uhat','U')
grid on
subplot(2,1,2)
%ax2 = nexttile;
%plot([1:length(yvals)],yvals,'r','lineWidth',1)
plot(t,yvals,'b','lineWidth',1)
xlabel ('Time')
ylabel('J')
grid on
%linkaxes([ax1 ax2],'x')
Telema Harry
Telema Harry el 27 de Nov. de 2020
I am trying to figure out why the results from my code is different.

Iniciar sesión para comentar.

Respuesta aceptada

VBBV
VBBV el 27 de Nov. de 2020
Editada: VBBV el 27 de Nov. de 2020
%true
xi = HPFnew.*sin((omega.*t+phase)*pi/180)
u = uhat +A*sin((omega*t+phase)*pi/180)
Try the above. Also see the filter function conditions

Más respuestas (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by