Using Help using filter function in Matlab

1 次查看（过去 30 天）
Telema Harry
Telema Harry 2020-11-27
编辑： Telema Harry 2020-11-29
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 个评论
显示 1更早的评论
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Digital Filter Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by