function phase_adj()
fsampl = 1000;      % частота дискретизации, кГц
freq = 68;          % Частота гармонического тона, кГц
len_buf = 4096;     % длина буфера, отсчёты
shift = pi/4;       % величина фазового сдвига между кусками, рад

% Заполняем буферы сигналом и аддитивным шумом
buffer_1 = sin(2*pi*(1:len_buf)*freq/fsampl)+(1e-5)*randn(1,len_buf);
buffer_2 = sin(2*pi*(1:len_buf)*freq/fsampl+shift)+(1e-5)*randn(1,len_buf);

% Пытаемся склеить без учёта фазовых соотношений
sum_buf_1 = [buffer_1, buffer_2];

% Делаем оценку спектров кусков и склейки, используя спектральные окна
specw_1=fft(blackman(len_buf)'.*buffer_1);
specw_2=fft(blackman(len_buf)'.*buffer_2);

sp_sum_1=fft(hamming(2*len_buf)'.*sum_buf_1);

% Находим фазовые спектры кусков
ang_1 = angle(specw_1);
ang_2 = angle(specw_2);

% Находим положение максимума модуля спектральной функции
[max_1, frq_1] = max(abs(specw_1));

% Находим фазу основного тона в конце 1-го куска (не доделано!)
ang_1_end = ang_1(frq_1);

% Находим разность фаз
diff_ang=ang_1_end-ang_2(frq_1);

% Находим поворачивающий множитель для частоты основного тона
rotator=-exp(i*diff_ang);

% Преобразуем второй буфер в комплексный вид
buf_2_hilb=hilbert(buffer_2);

% Разворачиваем его на величину сдвига фаз (нужно ещё ввести компенсацию по
% частоте, проще говоря - временнОй сдвиг)
buf_2_rot=real(buf_2_hilb*rotator);

% Склеиваем
sum_buf_2 = [buffer_1, buf_2_rot];

% Находим модуль спектра склейки
sp_sum_2=fft(hamming(2*len_buf)'.*sum_buf_2);

% Выводим результаты
figure (1)
plot (buffer_1(1:100))
hold on
plot (buffer_2(1:100), 'r')
hold off
grid on
title ('Signals')

figure (2)
plot ((4032:4160), sum_buf_1(4032:4160))
grid on
title ('Glue 1')

figure (3)
plot ((4032:4160), sum_buf_2(4032:4160), 'Color', [0.9 0 0])
grid on
title ('Glue 2')

figure (4)
semilogy (500:600, abs(sp_sum_1(500:600)))
hold on
semilogy (500:600, abs(sp_sum_2(500:600)), 'r')
hold off
grid on
title ('Spectra')


return;

function phase_adj()  
fsampl = 1000;      % частота дискретизации, кГц
freq = 40;          % частота гармонического тона, кГц

len_buf = 4096;     % длина буфера, отсчёты
shift = pi;       % величина фазового сдвига между кусками, рад

% Заполняем буферы сигналом и аддитивным шумом
buffer_1 = sin(2*pi*(0:len_buf-1)*freq/fsampl)+(1e-5)*randn(1,len_buf);
buffer_1 = buffer_1+0.03*sin(2*pi*(0:len_buf-1)*(freq+0.5)/fsampl);

shift_1=shift; % не корректно; токмо для ускорения процесса

buffer_2 = sin(2*pi*(0:len_buf-1)*freq/fsampl+shift)+(1e-5)*randn(1,len_buf);
buffer_2 = buffer_2+0.03*sin(2*pi*(0:len_buf-1)*(freq+0.5)/fsampl+shift_1);

% Пытаемся склеить без учёта фазовых соотношений
sum_buf_1 = [buffer_1, buffer_2];

% Делаем оценку спектров кусков и склейки, используя спектральные окна
specw_1=fft(hamming(len_buf)'.*buffer_1);
specw_2=fft(hamming(len_buf)'.*buffer_2);
sp_sum_1=fft(hamming(2*len_buf)'.*sum_buf_1);

% Находим фазовые спектры кусков
ang_1 = angle(specw_1);
ang_2 = angle(specw_2);

% Находим положение максимума модуля спектральной функции
[max_1, frq_1] = max(abs(specw_1));

% Находим разность фаз
diff=ang_1(frq_1)-ang_2(frq_1);
dif_ang=mod(diff+2*pi*freq/fsampl*len_buf, 2*pi);

% Находим поворачивающий множитель для частоты основного тона
rotator=exp(i*dif_ang);

% Преобразуем второй буфер в комплексный вид
buf_2_hilb=hilbert(buffer_2);

% Разворачиваем его на величину сдвига фаз (более правильно делать сдвиг во времени)
buf_2_rot=real(buf_2_hilb*rotator);

% Склеиваем
sum_buf_2 = [buffer_1, buf_2_rot];

% Находим модуль спектра склейки
sp_sum_2=fft(hamming(2*len_buf)'.*sum_buf_2);

% Выводим результаты
figure (1)
plot (1:100, buffer_1(1:100))
hold on
plot (1:100, buffer_2(1:100), 'r')
hold off
ylim([-1.5 1.5]);
grid on
title ('Signals')

figure (2)
plot ((4046:4146), sum_buf_1(4046:4146))
ylim([-1.5 1.5]);
grid on
title ('Glue 1')

figure (3)
plot ((4046:4146), sum_buf_2(4046:4146), 'Color', [0.9 0 0])
ylim([-1.5 1.5]);
grid on
title ('Glue 2')

figure (4)
semilogy (1:1000, abs(sp_sum_2(1:1000)))
grid on
title ('Spectrum of long sequence')

figure (5)
semilogy (1:1000, abs(specw_1(1:1000)))
grid on
title ('Spectrum of short sequence')
return;

function phase_adj()
fsampl = 1000;      % частота дискретизации, кГц
freq = 41;          % Частота гармонического тона, кГц
len_buf = 4096;     % длина буфера, отсчёты
gap = 2.5;          % интервал между кусками, мс, не более 8

% Формируем модель "реального" сигнала
signal = sin(2*pi*(0:4*len_buf)*freq/fsampl)+(1e-2)*randn(1, 4*len_buf+1)+...
     0.03*cos(2*pi*(0:4*len_buf)*(freq+0.5)/fsampl);

% Заполняем буферы (теперь всё точно)
buffer_1 = signal(1:len_buf);
sampl_gap = floor(1e6*gap/fsampl);
buffer_2 = signal(len_buf+1+sampl_gap:len_buf+sampl_gap+len_buf);

% Пытаемся склеить без учёта фазовых соотношений
sum_buf_1 = [buffer_1, buffer_2];

% Делаем оценку спектров кусков, используя спектральные окна
specw_1 = fft(chebwin(len_buf, 60)'.*buffer_1);
specw_2 = fft(chebwin(len_buf, 60)'.*buffer_2);

% Находим фазовые спектры кусков
ang_1 = angle(specw_1);
ang_2 = angle(specw_2);

% Находим положение максимума модуля спектральной функции
[max_1, frq_1] = max(abs(specw_1));

% Находим разность фаз с учётом набега за кусок
diff = ang_1(frq_1)-ang_2(frq_1);
dif_ang = mod(diff+2*pi*len_buf*freq/fsampl, 2*pi);

% Находим поворачивающий множитель для частоты основного тона
rotator = exp(i*dif_ang);

% Преобразуем второй буфер в комплексный вид
buf_2_hilb=hilbert(buffer_2);

% Разворачиваем его на величину сдвига фаз (более правильно делать сдвиг во времени)
buf_2_rot=real(buf_2_hilb*rotator);

% Склеиваем
sum_buf_2 = [buffer_1, buf_2_rot];

% Находим модуль спектра склейки
sp_sum_2=fft(chebwin(2*len_buf, 60)'.*sum_buf_2);

% Выводим результаты
bins = (0:len_buf-1)*fsampl/len_buf;
figure (1)
plot (1:100, buffer_1(1:100))
hold on
plot (1:100, buffer_2(1:100), 'r')
hold off
ylim([-1.5 1.5]);
grid on
title ('Signals')

figure (2)
plot ((4046:4146), sum_buf_1(4046:4146))
ylim([-1.5 1.5]);
grid on
title ('Glue 1')

figure (3)
plot ((4046:4146), sum_buf_2(4046:4146), 'Color', [0.9 0 0])
ylim([-1.5 1.5]);
grid on
title ('Glue 2')

figure (4)
semilogy (bins(1:1000), abs(specw_1(1:1000)))
ylim([0.1 1000]);
grid on
title('Spectrum of short sequence')

figure (5)
semilogy (bins(1:1000)/2, abs(sp_sum_2(1:1000)))
ylim([0.1 1000]);
grid on
title('Spectrum of long sequence')
return;