陷波器程序如下，麻烦帮忙看看

clc;

clear all;

close all;

x = 0 : 0.05 : 10;

s= 210 * exp(-0.25 * x) .* cos(2 * pi * 1.5 * x + pi) + 110 * exp(0.25 * x) .* cos(2 * pi * 1.0 * x + pi/2 );

%load data;

%s=cc;

N=length(s);                            % 信号长度

fs=1000;                                % 采样频率

n=1:N;

n2=1:N/2;

tt=(n-1)/fs;                            % 时间刻度

ff=(n2-1)*fs/N;                         % 频率刻度

X=fft(s);                               % 谱分析

figure(1)

subplot 311;

plot(tt,s,'k');

title('原始数据'); xlabel('时间/s'); ylabel('幅值');

X=X/max(abs(X));

subplot 312;

plot(ff,abs(X(n2)),'k');

axis tight; title('谱分析');

xlabel('频率/Hz'); ylabel('幅值');

subplot 313;

plot(tt,abs(X(n)));

% for k=1:5                             % 自适应陷波器

%     j=(k-1)*2+1;                        % 设置50Hz和它的奇次谐波频率

%     f0=116.73*j;

%     j=(k-1)*2+1;                        % 设置50Hz和它的奇次谐波频率

f0=116.73;

x1=cos(2*pi*tt*f0);                 % 设置x1和x2

x2=sin(2*pi*tt*f0);

w1=0;                               % %初始化w1和w2

w2=100;

e=zeros(1,N);                       % %初始化e和y

y=zeros(1,N);

mu=0.0000001;                             % 设置迭代步长

for i=1:N                           % 自适应陷波器

y(i)=w1*x1(i)+w2*x2(i);         % 计算y

e(i)=s(i)-y(i);                 % 计算e

w1=w1+mu*e(i)*x1(i);            % 调整w

w2=w2+mu*e(i)*x2(i);

end

x=e;

% end

output1=e;                               % 陷波器输出