Using Help using filter function in Matlab

2 次查看(过去 30 天)
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 个评论
Telema Harry
Telema Harry 2020-11-27
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 2020-11-27
I am trying to figure out why the results from my code is different.

请先登录,再进行评论。

采纳的回答

VBBV
VBBV 2020-11-27
编辑:VBBV 2020-11-27
%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

更多回答(0 个)

Community Treasure Hunt

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

Start Hunting!

Translated by