Гетеродинное преобразование квадратурного сигнала и доплеровское смещение

% clc;
fs = 102e6; % частота дискретизации
fh = 1530e6; %частота гетеродинного преобразования
fc = 1602e6; % несущая частота ГЛОНАСС
collect_number = 20; %число накоплений
N = collect_number*fs/1000;
fileID = fopen('D:\no_jamm_2_10_10_0_0.txt');
C = textscan(fileID,'%f %f %f %f %f %f %f %f', N);
data = [C{1,:}];
fclose(fileID);
ch_adc=[ 1 ]; %берем первый вход АЦП
 
I = zeros(1,max(size(data))); %формируем массив размером 1 х максимальную длину данных
Q = zeros(1,max(size(data)));
for i = ch_adc
    I = I + (data(:,2*i-1)-mean(data(:,2*i-1)))';
    Q = Q + (data(:,2*i)-mean(data(:,2*i)))';
end
IQ = I + 1i*Q;
 
fM = 0.511*10^6; %скорость передачи данных
fs_fM=fs/fM; %отношение частоты дискретизации к символьной скорости
C = 7.9e3; %скорость спутника в м/с(первая космическая)
t = 1/fs; %дискретное время
fi = fc - fh; %вычисляем промежуточную частоту(замена гетеродинному преобразованию)
 
S = I.*cos(2*pi*fi*t) + i*Q.*sin(2*pi*fi*t); %формируем сигнал
 
 
S_h = S.*sin(2*pi*fi*t); %гетеродинное преобразование

fi = fc - fh; %вычисляем промежуточную частоту
 
S = IQ.*(sin(2*pi*fi*t)+cos(2*pi*fi*t)); %формируем сигнал
 
 
S_h = S.*exp(2*i*pi*fi*t); %гетеродинное преобразование

S_h = S.*exp(2*i*pi*fi*t); %гетеродинное преобразование

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% ПРИБЛИЗИТЕЛЬНЫЙ ПОИСК ДОПЛЕРОВСКОГО СДВИГА
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for k=0:100:400
    F = fi + k; % приблизительный поиск частоты доплера
    Z = Y.*cos(2*pi*(fi+k)*t); %гетеродинное преобразование с промежуточной + поиск доплеровской частоты
    
    Num = [1,2,1]; %коэффициенты числителя 1-ой секции ФНЧ
    Den = [1,-1.9485228131107961,0.95098699434575018]; %коэффициенты знаменателя 1-ой секции
    y = filter(Num, Den, Z); %фильтрация по 1-ой секции
    
    Num = [1,1,0]; %коэффициенты числителя 2-ой секции ФНЧ
    Den = [1,-0.9509567815010681,0]; %коэффициенты знаменателя 2-ой секции
    X = filter(Num, Den, y); %фильтрация по 2-ой секции
    
    F1 = fft(M(fix(fM*t)+1));
    Kf = conj(F1);
    F2 = fft(X);
    Fv = F2.*Kf;
    Sv = ifft(Fv)
    
    figure
    subplot(2,1,1)
    plot(Sv)
    grid on
    
    subplot(2,1,2)
    plot(X)
    grid on
end