Ошибка в матлабе!

% Исходные данные
U_m=input('Амплитуда сигнала (В): ');          % Амплитуда сигнала
F_prom=63;                                     % Промежуточная частота (МГц)
T_rep=10;                                      % Период повторения импульсов (мкс)
F_rep=1/T_rep;                                 % Частота повторения (МГц)
tau_imp=1;                                     % Длительность импульса (мкс)
N_imp=input('Количество импульсов в пачке: '); % Количество импульсов в пачке
T_pack=N_imp*T_rep;                            % Длительность пачки радиоимпульсов (мкс)
F_Doppler=input('Доплеровский сдвиг частоты (кГц): ')*10^(-3); % Доплеровский сдвиг частоты (МГц)
fs=84;
dt=1/fs;                                       % Шаг дискретизации (мкс)
tay_delay=input('Задержка отраженного сигнала (мкс): '); % Задержка отраженного сигнала (мкс)
N_sample_1=T_rep*fs;                           % Количество отсчетов в одном периоде
N_sample_imp=tau_imp*fs;                       % Количество отсчетов в импульсе
 
figure(1)
% Отр. от цели сигнал. Сигнал представляет собой ког. пачку импульсов.
t=0:0.001:T_pack;                              % Интервал времени пачки импульсов
d=[0.5:10:T_pack];                             % Массив значений центров импульсов
y=U_m*pulstran(t-tay_delay,d,'rectpuls',1);    % Формирование последовательности прямоугольных импульсов
b=cos(2*pi*(F_prom+F_Doppler)*t);              % Формирование гармонического сигнала т.е. ВЧ заполнение в импульсе
x=y.*b;                                        % Получение радиоимпульса путем умножения огибающей сигнала и гармонического сигнала
plot(t,x),grid                                 % График принятого отраженного сигнала
ylim([-U_m-0.5,U_m+0.5])                       % Интервал значений по оси y
 
figure(2)
% Аналоговый полосовой фильтр f0=63 МГц, полоса пропускания 14 МГц.
Center_frequency=63;                           % Центральная частота аналогового полосового фильтра
[b,a]=butter(5,[2*pi*(Center_frequency-7) 2*pi*(Center_frequency+7)],'s'); % Синтез коэффициентов фильтра
f=0:0.1:90;                                    % Массив значений частоты (МГц)
h=freqs(b,a,2*pi*f);                           % Комплексный коэффициент передачи
plot(f,abs(h))                                 % График АЧХ полосового фильтра                                 
xlim([30 90])                                  % Интервал значений по оси x
H=tf(b,a);                                     % Передаточная функция соответствующая своим коэффициентам числителя и знаменателя
 
figure(3)
% Формирование шума. Нормальный случайный процесс.
T=0:dt:T_pack-dt;                              % Массив дискретного времени (мкс)
Shum=randn(1,length(T));                       % Формирование шума
Shuml=lsim(H,Shum,T);                          % Шум на выходе аналогового фильтра
Shum2=Shuml';                                  % Изменение вектора столбца на вектор строку
plot(T,Shum2),grid                             % График шума на выходе аналогового фильтра
 
figure(4)
% Оцифрованный (отраженный) сигнал. Сигнал на выходе АЦП
y_digitize=U_m*pulstran(T-tay_delay,d,'rectpuls',1);   % Формирование последовательности прямоугольных импульсов
b_digitize=cos(2*pi*(F_prom+F_Doppler)*T);             % Формирование гармонического сигнала т.е. ВЧ заполнение в импульсе
x_digitize=y_digitize.*b_digitize+Shum2;               
stem(T,x_digitize),grid 
 
figure(5)
% Вычисление спектра от дискретного (отраженного) сигнала.
NFFT=N_sample_1*N_imp;                         % Количество дискретных отсчетов сигнала
df=fs/NFFT;                                    % Шаг по частоте в спектре (МГц)
f=(0:NFFT-1)*df;                               % Диапозон частот в спектре (МГц)
Spectr=abs(fft(x_digitize,NFFT));              % Вычисление спектра
plot(f,Spectr),grid                            % График дискретного спектра
xlim([0,fs])                                   % Интервал значений по оси x

T=0:dt/10:T_pack-dt;