Быстрое преобразование Фурье

              
%% начальные данные и сигнал
clear;
t=0:0.000008:0.0004; %0.000000001
sigma_1=0;
sigma_2=0;
omega_1=50000;
omega_2=50000;    %400мкс   250мкс  53мкс-задержка
 
 
tz=0.000053;
a1=sin(2*pi*omega_2*t);
a2=sin(2*pi*omega_2*(t-tz)); 
y1=a1.*(t>0.00002 & t<0.00027)+sigma_2*rand(size(t));
y2=a2.*(t>0.00002+tz & t<0.00027+tz)+sigma_2*rand(size(t));
 
Fd=125000;% Частота дискретизации (Гц)
Epsilon=8000;% Область максимума
Ff=64;% Количество линий Фурье спектра  65536
Fe=Fd/Ff;  
 
subplot(5,1,1);
plot(y1);
subplot(5,1,2);
plot(y2);
 
%% Построение и анимация взаимокорреляции
C=0; j=1; k=1;
 
y3=xcorr(y1,y2);
  [C1,j1] = max(y3);
 
subplot(5,1,3);
plot(y3);
 
 F=0:Fd/Ff:Fd/2-1/Ff;% Массив частот вычисляемого спектра Фурье
 
FftS1=fft(y1,Ff);% Амплитуды преобразования Фурье сигнала
FftS2=fft(y2,Ff);% Амплитуды преобразования Фурье сигнала
FftS=FftS1.*conj(FftS2);
FftSc=2*abs(FftS)./Ff;
FftSh=unwrap(angle(FftS));
[C,i] = max(FftSc);
i=ceil(50000/Fe);
 
subplot(5,1,4);
plot(F(max(1,i-Epsilon):min(length(F),i+Epsilon)),...
        FftSc(max(1,i-Epsilon):min(length(F),i+Epsilon)),'-r',...
        F(i),FftSc(i),'ko');
subplot(5,1,5);
plot(F(max(1,i-Epsilon):min(length(F),i+Epsilon)),...
        FftSh(max(1,i-Epsilon):min(length(F),i+Epsilon)),'-r',...
        F(i),FftSh(i),'ko');
 
  
a1=diff(FftSh((i-1):(i+1)));
mean(a1)/(2*pi*Fe)
 
figure;
subplot(2,1,1);
plot(y1);
subplot(2,1,2);
plot(y2);
 
%% БПФ ФВК
Epsilon=100; Ff=128;
F=0:Fd/Ff:Fd/2-1/Ff;% Массив частот вычисляемого спектра Фурье
FfVk1=fft(y3,Ff);
FfVk=conj(FfVk1);
FfSc=2*abs(FfVk)./Ff;
FfSh=unwrap(angle(FfVk));
Fe=Fd/Ff; 
i=ceil(50000/Fe);
 
figure;
subplot(2,1,1);
plot(F(max(1,i-Epsilon):min(length(F),i+Epsilon)),...
        FfSc(max(1,i-Epsilon):min(length(F),i+Epsilon)));
subplot(2,1,2);
plot(F(max(1,i-Epsilon):min(length(F),i+Epsilon)),...
        FfSh(max(1,i-Epsilon):min(length(F),i+Epsilon)));
 
a2=diff(FfSh((i-2):(i+1)));
mean(a2)/(2*pi*Fe)