spectrumanalysis.m

本人编写的信号滤波并重构

function y=positionalg()
load data3.mat -mat;
mylength=3e3;
pos=[];
pos1=[];
pos2=[];

% subplot(2,1,1),plot(tsignal);
% subplot(2,1,2),plot(data2(1:mylength,1),'r');
for i=55:70
s1=data3((i-1)*mylength+1:i*mylength,1);
s2=data3((i-1)*mylength+1:i*mylength,2);
s1=(s1-mean(s1))/sqrt(sum(s1.^2)/mylength);
s2=(s2-mean(s2))/sqrt(sum(s2.^2)/mylength);
% s1=smooth(s1,5,'rlowess');
% s2=smooth(s2,5,'rlowess');
hold on;
plot(s1);
plot(-s2,'r');
s1=[zeros(mylength/2,1);s1;zeros(mylength/2,1)];
s2=[zeros(mylength/2,1);s2;zeros(mylength/2,1)];
fsignal1=fft(s1);
fsignal2=fft(s2);
% for i=1:2*mylength
%     if abs(fsignal1(i))<3
%         fsignal1(i)=0;
%     end
%     if abs(fsignal2(i))<3
%         fsignal2(i)=0;
%     end
% end
% fsignal1=[zeros(mylength/2,1);fsignal1;zeros(mylength/2,1)];
% fsignal2=[zeros(mylength/2,1);fsignal2;zeros(mylength/2,1)];
fsignal=fsignal1.*conj(fsignal2);
tsignal=abs(fftshift(ifft(fsignal)));
subplot(2,2,1),plot(tsignal);
[out,k]=max(tsignal);
pos=[pos (mylength-k)];
% s1=ifft(fsignal1);
% s2=ifft(fsignal2);
% 
% [A,d1,Ea,Ey,P] = kf(d11,d11,2,1);
% [A,d2,Ea,Ey,P] = kf(d21,d21,2,1);

% Fs = 1000000;  % Sampling Frequency
% 
% Fstop = 10;           % Stopband Frequency
% Fpass = 50;           % Passband Frequency
% Dstop = 0.01;           % Stopband Attenuation
% Dpass = 0.057501127785;  % Passband Ripple
% flag  = 'noscale';       % Sampling Flag
% 
% % Calculate the order from the parameters using KAISERORD.
% [N,Wn,BETA,TYPE] = kaiserord([Fstop Fpass]/(Fs/2), [0 1], [Dpass Dstop]);
% 
% % Calculate the coefficients using the FIR1 function.
% b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
% Hd = dfilt.dffir(b);
% 
% 
% d1=filter(Hd,s1);
% d2=filter(Hd,s2);


y=abs(xcov(s1,s2));
[out,k]=max(y);
pos2=[pos2 (2*mylength-k)];
hold on
subplot(2,2,3), plot(y,'-r');
hold on
end
%hold off
subplot(2,2,4),plot(ones(size(pos)),pos,'o',2*ones(size(pos1)),pos1,'o',3*ones(size(pos2)),pos2,'o');