clc
clear all
%% FFT采样
fs=3200;
N=128; %采样频率和数据点数
n=0:N-1;
t=n/fs; %时间序列
%x=98.23*sin(2*pi*50*t+(1/6)*pi)+81*sin(2.1*2*pi*50*t+(1/3)*pi)+60*sin(5*2*pi*50*t-(1/6)*pi)+35*sin(7.3*2*pi*50*t-(1/3)*pi)+20*sin(9*2*pi*50*t-(1/4)*pi)+9*sin(11.5*2*pi*50*t+(1/4)*pi); %信号
x=0.8*sin(4.5*2*pi*50*t+pi/4)+0.3*sin(2*2*pi*50*t+pi/5); %信号
y=fft(x,N); %对信号进行快速Fourier变换
f=fs*(0:length(y)/2)/length(y);
z=angle(y)*180/pi;
mag=abs(y);
fz=mag(1:length(y)/2+1)*2/N; %求得Fourier变换后的振幅
figure(1)
plot(f,fz); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');
title('FFT采样');
grid on;
%% BP神经网络初始化
inputNums=1;
outputNums=1;
hideNums=4;
maxcount=3000;
sampleNums=N;
precision=1e-9;
xxl1=0.02;
xxl2=250;
count=1;
error=zeros(1,count);
errorp=zeros(1,sampleNums);
w1=rand(inputNums,hideNums/2);
w2=rand(inputNums,hideNums/2);
b=[1*pi*200,1*pi*120];
%% BP神经网络训练
while (count<=maxcount)
c=1;
while (c<=sampleNums)
for j=1:hideNums/2
net1(j)=cos(b(j)*t(c));
net2(j)=sin(b(j)*t(c));
end
for j=1:hideNums/2
out(j)=net1(j)*w1(j)+net2(j)*w2(j);
end
sc(c)=sum(out);
errortmp=0.0;
errortmp=errortmp+(x(c)-sc(c))^2;
e(c)=x(c)-sc(c);
errorp(c)=0.5*errortmp;
for i=1:hideNums/2
dw1(i)=xxl1*e(c)*net1(i);
dw2(i)=xxl1*e(c)*net2(i);
db(i)=xxl2*e(c)*t(c)*(-w1(i)*net2(i)+w2(i)*net1(i));
end
for k=1:hideNums/2
w1(k)=w1(k)+dw1(k);
w2(k)=w2(k)+dw2(k);
b(k)=b(k)+db(k);
end
c=c+1;
end
double tmp;
tmp=0.0;
for i=1:sampleNums
tmp=tmp+errorp(i)*errorp(i);
end
tmp=tmp/c;
error(count)=sqrt(tmp);
if (error(count)<precision)
break;
end
count=count+1;
end
%% 误差精度曲线
figure(2)
p=1:count;
plot(p,error(p),'-')
xlabel('迭代次数');
ylabel('误差');
title('BP误差曲线');
%% 谐波分析结果
for i=1:hideNums/2
Amp(i)=sqrt(w1(i)^2+w2(i)^2);
Phase(i)=atan(w1(i)/w2(i))/pi*180;
HFre(i)=b(i)/2/pi;
end
Amp
Phase
HFre
%%
