Форум кафедры Техники и Электрофизики Высоких Напряжений

ti = 5e-5; % параметры импульса
tf = 1e-5;
l = 45000; % Длина линии, м
tau = l/3e8; % время пробега по линии без потерь со скоростью света
N=2^22; % число точек дискретизации
fkr=2e7; % критическая частота, Гц

dt = 1/(2*fkr); % шаг дискретизации
t = 0:dt:dt*(N-1); % вектор-строка отсчетов по времени, с
Uvh = exp(-t/ti)-exp(-t/tf); % исходный сигнал, В

Tk = dt*(N-1); % область решения (размер временного окна), с
df = 1/Tk; % шаг по частоте, Гц
n = -N/2+1:1:N/2; % индексы отсчетов по частоте
f = n*df; % вектор-строка отсетов по частоте
UVH = fftshift(fft(Uvh))*dt; % спектр исходного сигнала, В*с

fr = [logspace(1,7,71) 2e7]; % набор частот для расчета коэффициентов распространения

% выделяем память под массивы коэффициентов распространения
Gkt0=zeros(1,length(fr));
Gkt1=zeros(1,length(fr));
txt=''; % строка индикатора прогресса расчета
for I=1:length(fr)
    % ----- вывод текущего прогресса расчета -----
    for jj=1:length(txt); % стираем символы, выведенные при предыдущем проходе цикла
        fprintf('\b');
    end
    clear txt;
    txt = ['Выполнено ' sprintf('%3.1f',I/length(fr)*100) '%'];
    fprintf('%s',txt);
    % ----- конец вывода текущего прогресса расчета -----
   
    [~, ~, Gkt, ~, ~, ~, ~]=p12(Data,fr(I),roz);
    Gkt1(I)=Gkt(1);
    Gkt0(I)=Gkt(3);
end
fprintf('\n');  % переводим курсор на новую строку

ffull = [-fliplr(fr) fr]; % расширяем набор частот в отрицательную область

% перекидываем коэффициент распространения нулевого канала в область отрицательных частот
Gfullr = [fliplr(real(Gkt0)) real(Gkt0)];
Gfulli = [-fliplr(imag(Gkt0)) imag(Gkt0)];
Gkfull = Gfullr+1j*Gfulli; % объединяем действительную и мнимую части
G0 = spline(ffull,Gkfull,f); % интерполируем на набор частот для ДПФ

% аналогично для воздушного канала
Gfullr = [fliplr(real(Gkt1)) real(Gkt1)];
Gfulli = [-fliplr(imag(Gkt1)) imag(Gkt1)];
Gkfull = Gfullr+1j*Gfulli;
G1 = spline(ffull,Gkfull,f);

K0 = exp(-G0*l); % передаточная функция для нулевого канала
K1 = exp(-G1*l); % передаточная функция для воздушного канала
Klossless = exp(-1j*2*pi*f*tau); % передаточная функция для линии потерь
% находим спектры импульса после пробега
UOUT0 = UVH.*K0;
UOUT1 = UVH.*K1;
UOUTlossless = UVH.*Klossless;
% делаем обратное дискретное преобразование Фурье
Uout0 = ifft(ifftshift(UOUT0))*N/Tk;
Uout1 = ifft(ifftshift(UOUT1))*N/Tk;
Uoutlossless = ifft(ifftshift(UOUTlossless))*N/Tk;
% строим графики
plot(t*1e6,([Uvh; Uout0; Uout1; Uoutlossless]));
xlim([0 500]);
ylim([0 0.6]);
xlabel('t, мкс')
ylabel('U, В')
legend('исходный импульс', 'пробег 45 км по нулевому каналу', 'пробег 45 км по воздушному каналу', 'пробег 45 км по линии без потерь')
grid on