Using Help using filter function in Matlab

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.

This is the result expected. I don't know what i'm doing wrong?

[Mathieu NOE]

hello

so what is the code that genrated the "right" plot and why yours differs ?

[Telema Harry]

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]

I am trying to figure out why the results from my code is different.

Answer 1

[VBBV]

%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

[Telema Harry]

This time, I got a different response.