Перейти к содержанию
    

Не стабилизируется ФАПЧ в Матлаб

Помогите пожалуйста найти ошибку в Матлаб коде петли ФАПЧ

%% Тест ФАПЧ 

clear all; 
close all; 

%% Исходные установки

%bit_stream = (rand(1, 20) > 0.5);
bit_stream = [0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
N = length(bit_stream);

% Частота семплирования 288 кГц (8 выборок на период)
fs = 288000;
% Несущая частота относительно частоты выборок 36 кГц
f = 0.125;
% Частота доплера 
fdop = 0;
% время выборки
Ts = 1/fs;
fprintf('Время выборки Ts %d\n', Ts);
% Амплитуда
A = 1.0;
% Символьная частота Гц
fsim = 3600;
% Количество выборок на период несущей 
ns = 1/f;
% Количество выборок необх. для одного символа 
% при символьной скорости fsim
tsymbol = fs*f*ns/fsim;
fprintf('Кол. пер. несущ. в одном симв. tsymbol %d\n', fs*f/fsim);
% несущая частота Гц
fref = f*fs;
fprintf('Несущая %d\n', fref);
ttotal = N*tsymbol;

%% Входной сигнал

t = 0;
for i = 1:length(bit_stream)
%  phase = pi * bit_stream(i);
 phase = 0;
 for j = 1:tsymbol
   % входной сигнал
   BPSK_signal(t+1) = cos(2*pi*(f+fdop)*t + phase); 
   t = t + 1;
 end 
end

%% Расчет коэффициентов петлевого фильтра

% натуральная частота Гц
fn = 500; 
% коэффициент усиления NCO (по факту минимальный шаг соотв. разр. сетке)  
Knco= 1/4096; 
% коэфф. усил. фазового детектора 1/cycles
KP = 2;              
% демпинг фактор 
zeta = 1.0;       
% круговая натуральная частота rad/s fn = fn/(2*pi)
wn = 2*pi*fn;      
% пропорциональный коэффициент петлевого фильтра
KL= 2*zeta*wn*Ts/(KP*Knco);      
% интегральный коэффициент петлевого фильтра
KI= wn^2*Ts^2/(KP*Knco);
% распечатка
fprintf('Пропорц. коэфф KL %d\n', KL);
fprintf('Интегр. коэфф KI %d\n', KI);

% Расчет коэффициентов передаточной функции с замкнутого контура u / ref_phase
% CL(z) = (b0 + b1z^-1)/(a2Z^-2 + a1z^-1 + 1)
b0= KP*KL*Knco;
b1= KP*Knco*(KI - KL);
a1= KP*KL*Knco - 2;
a2= 1 + KP*Knco*(KI - KL);
% коэффициенты петлевого фильтра (numerator denominator)  
b = [b0 b1];        
a = [1 a1 a2];       
fprintf('b %d\n', B);
fprintf('a %d\n', a);

%% Собственно ФАПЧ

% начальная частота NCO Гц 
fnco = fref; 
fprintf('fnco Гц %d\n', fnco);
% индексы времени модели
n = 0:ttotal; 
% циклы начальной фазы опорного сигнала 
init_phase = 0.7; 
% циклы фазы опорного сигнала
ref_phase = fref*n*Ts + init_phase; 
% циклы фазы опорн. сигн. по mod 1
ref_phase = mod(ref_phase,1); 
fprintf('ref_phase %d\n', ref_phase);

u(1) = 0;
ur(1) = 0;
int(1)= 0;
% начальная фазовая ошибка
phase_error(1) = -init_phase;
% начальная значение на входе NCO
vtune(1) = -init_phase*KL;
fprintf('vtune(1) %d\n', vtune(1));

%per = 0.05;

for step = 2:ttotal

 % NCO

 x = fnco*Ts + u(step-1) + vtune(step-1)*Knco;  % циклы NCO фазы
%  x = fnco*Ts + u(step-1);       % циклы NCO фазы
 u(step) = mod(x,1);             % циклы NCO фазы по mod 1
 s = cos(2*pi*u(step-1));        % NCO sin выход
 y(step)= round(2^15*s)/2^15;    % квантованный выход синуса NCO

 % Fref входной сигнал

%  xr = fref*Ts + ur(step-1);       % циклы fref
%  ur(step) = mod(xr,1);            % циклы fref фазы по mod 1
%  sr = sin(2*pi*ur(n-1));         % sin выход fref 
 sr = BPSK_signal(step-1);        % sin выход fref 
 yr(step)= round(2^15*sr)/2^15;   % квантованный выход fref

 % выход умножителя опорного и NCO
 Detect(step) = yr(step)*y(step);

 % фазовый детектор

 per= ref_phase(step-1) - u(step-1); % ошибка фазы
 per= 2*(mod(per+1/2,1) - 1/2);      % обертка, если пересечение фаз +/- 1/2 цикла
 phase_error(step) = per;

 % Петлевой фильтр
 int(step) = KI*per + int(step-1);   % интегратор
 vtune(step) = int(step) + KL*per;   % выход петлевого фильтра

end  

%% графика

% фазовая ошибка с фазового детектора
figure
plot(phase_error, 'b-', 'LineWidth', 2),grid
axis([0 3 -1 1])
xlabel('t (ms)'),ylabel('phase_error'),title('фазовая ошибка с фаз. детектора')

figure;
plot(Detect, 'r-', 'LineWidth', 2);
hold on;
grid on;
title('Выход умножителя входного и NCO сигналов')

% выход VCO
figure
plot(vtune, 'r-', 'LineWidth', 2),grid
%axis([0 30 -1 5])
title('выходной сигнал петлевого фильтра')

% sin NCO и fref
figure;
plot(yr, 'b-', 'LineWidth', 2);
hold on;
grid on;
plot(y, 'r-', 'LineWidth', 2);
hold on
hold off;
title('Входной и опорный сигналы')

% Частотная характеристика замкнутого контура
figure
u = 0:.1:.9;
f= 10* 10 .^u;                  % log-scale frequencies
f = [f 10*f 100*f 1000*f];
z = exp(j*2*pi*f/fs);           % complex frequency z
CL= (b0 + b1*z.^-1)./(1 + a1*z.^-1 + a2*z.^-2);    % closed-loop response
CL_dB= 20*log10(abs(CL));
semilogx(f,CL_dB),grid
xlabel('Hz'),ylabel('CL(z) dB'),title('Частотная характеристика замкнутого контура')

Фапч посчитана по правилам приведенным в букваре pll_notes

Посчитаны интегрирующий и пропорциональный коэффициенты.

Сделана собственно петля с фазовым детектором, интегратором, и NCO.

Но, блин не выходит частота NCO на частоту входного сигнала, даже при нулевой расстройке. Не врубаюсь почему.

Для контроля просто вывел на картинку выход с умножителя опорного на NCO

Вот что получается:

post-39850-1524748249_thumb.jpg

Изменено пользователем Acvarif

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Нарисуйте уже в симулинке вашу петлю из простейших блоков.

Пока смутно представляю как это сделать, поскольку в симулинк сильно плаваю.

Но ошибку нашел. ФАПЧ работает, неплохо работает. Таки да, ей без разницы какой сигнал, модулированный или нет.

На картинках спектры выхода NCO и входного сигнала с шумом с соотн 1, при долере 120 Гц на несущей 36 000 Гц.

post-39850-1524824469_thumb.jpg

Остаются еще неясности.

Например:

1. Как можно просто пояснить что такое натуральная частота?

2. Нет полного представления как петлю пристроить к BPSK демодулятору.

Изменено пользователем Acvarif

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Плавайте дальше. Вам даже картинку нарисовали в книжке. Что смутного в сумматоре, умножителе, задержке? Если с 2x2 проблема, стоит ли вообще браться за такие задачи?

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Плавайте дальше. Вам даже картинку нарисовали в книжке. Что смутного в сумматоре, умножителе, задержке? Если с 2x2 проблема, стоит ли вообще браться за такие задачи?

Не в этом дело. Любую среду нужно освоить, хотя-бы механически, сначала на простейших примерах и т. д. Что я и делаю...

Вы же не будете спорить, что зная как провести один проводник к другому нельзя сразу развести плату в PCAD.

Или нарисовав схему запустить ее моделирование в ModelSim.

Даже неплохо владея С++ нельзя сразу писать код в Матлаб.

Изменено пользователем Acvarif

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

ФАПЧ работает, неплохо работает. Таки да, ей без разницы какой сигнал, модулированный или нет.

На картинках спектры выхода NCO и входного сигнала с шумом с соотн 1, при долере 120 Гц на несущей 36 000 Гц.

 

Acvarif, Вы несколько раз приводили предельные значения SNR < 1 входного сигнала, при которых работает Ваша система. Можете привести подробный расчет ( сигнал+шум) с учетом полосы занимаемой сигналом, который использовали в моделировании.

Несущая Fsig=36000 Гц; Амплитуда Asig=1 (???).

Символьная Fsim=Fsig/n=100 Гц (???). (n=360 ???).

Частота дискретизации Fsr=288000 Гц. ( ???)

Эффективная полоса сигнала dFsig=Fsim (2Fsim ... другая ???) Гц.

Формла расчета ( сигнал+шум) = .....???

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Acvarif, Вы несколько раз приводили предельные значения SNR < 1 входного сигнала, при которых работает Ваша система. Можете привести подробный расчет ( сигнал+шум) с учетом полосы занимаемой сигналом, который использовали в моделировании.

Несущая Fsig=36000 Гц; Амплитуда Asig=1 (???).

Символьная Fsim=Fsig/n=100 Гц (???). (n=360 ???).

Частота дискретизации Fsr=288000 Гц. ( ???)

Эффективная полоса сигнала dFsig=Fsim (2Fsim ... другая ???) Гц.

Формла расчета ( сигнал+шум) = .....???

Нет. С учетом полосы занимаемой сигналом нет. Полоса BPSK сигнала равна удвоенной символьной частоте, тоесть 200 Гц

В модели просто смешиваю сигнал с шумом 1:1 с амплитудой сигнала 1. Очевидно что это не совсем верно.

Буду признателен если подскажете как правильно сделать расчет.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Нет. С учетом полосы занимаемой сигналом нет. Полоса BPSK сигнала равна удвоенной символьной частоте, тоесть 200 Гц

В модели просто смешиваю сигнал с шумом 1:1 с амплитудой сигнала 1. Очевидно что это не совсем верно.

Буду признателен если подскажете как правильно сделать расчет.

 

Эффективная полоса сигнала в 200 Гц это по первым нулям АЧХ. Обычно используют полосу содержащуюю половину мощности. Но это уже детали.

1. Пропустить сформированный сигнал BPSK через полосовой фильтр dFsig, и вычислить СКЗ напряжения вашего сигнала (Usig). Можно вычислить через fft, используя только нужные частоты попадающие в полосу dFsig, Или рассчитать это значение на счетах.

2. Сформировать шум snoise=randn(size(signal)); имеющий единичный уровень СКЗ в полосе (Fsr/2).

3. Рассчитать уровень шума попадающий в полосу вашего сигнала

UnoiseFilt=sqrt(2*dFsig/Fsr);

4. Вычислить SignalNoise=signal+snoise.*(Usig/(UnoiseFilt*SNR));

 

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Эффективная полоса сигнала в 200 Гц это по первым нулям АЧХ. Обычно используют полосу содержащуюю половину мощности. Но это уже детали.

1. Пропустить сформированный сигнал BPSK через полосовой фильтр dFsig, и вычислить СКЗ напряжения вашего сигнала (Usig). Можно вычислить через fft, используя только нужные частоты попадающие в полосу dFsig, Или рассчитать это значение на счетах.

2. Сформировать шум snoise=randn(size(signal)); имеющий единичный уровень СКЗ в полосе (Fsr/2).

3. Рассчитать уровень шума попадающий в полосу вашего сигнала

UnoiseFilt=sqrt(2*dFsig/Fsr);

4. Вычислить SignalNoise=signal+snoise.*(Usig/(UnoiseFilt*SNR));

Спасибо. В общих чертах вроде понятно.

Поскольку сталкиваюсь с этим впервые то на деле, сходу не получится.

Что такое полоса содержащая половину мощности?

1. Предполагается что приемный усилитель будет имень в своем составе полосовой фильтр с центр. частотой 36 кГц. Вопрос в том

какая у него должна быть полоса. Если заложить всего 200 Гц это будет неверно, - высокая селективность это плохая устойчивость для фазоманипулированного сигнала. Если заложить всю возможную полосу 2000 Гц которую обеспечивает гидроакустический преобразователь тоже будет не верно. Лишние спектральные компоненты в приемн. тракте тоже не нужны.

СКЗ напряжения вашего сигнала - нет представления что это ...

2. Матлаб примеры предлагают 0.00001*randn(1,N); где N кол. выборок

Что такое полоса (Fsr/2)?

3. Вроде понятно

4. Тоже понятно.

Матлаб предлагает вычислять SNR так https://www.mathworks.com/help/signal/ref/snr.html

где все завязано на частоте сигнала, кол. выборок на период, длительности сигнала.

 

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Эффективная полоса...

 

Вот смотрю графики BER в книжках и не вижу там никакой полосы.

 

1.Нужно измерить мощность сигнала.

2. Задать параметры блока AWGN Channel https://www.mathworks.com/help/comm/ref/awgnchannel.html

Нужный Eb/No.

Количество бит передаваемых на символ.

Мощность сигнала измеренную.

Длительность символа.

3. Практически измерить BER при заданном Eb/No и убедиться, что совпадает с теорией.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Вот смотрю графики BER в книжках и не вижу там никакой полосы.

BER это что? Битовая ошибка?

Полоса BPSK http://www.dsplib.ru/content/bpsk/bpsk.html

В матлаб SNR оцениваю так

% отношение сигнал/шум в децибелах
SNR = 0;
% смесь мод.сигнала с шумом
BPSK_signal_awgn = awgn( BPSK_signal_t, SNR);

% фактическая мощность сигнала 
Sbpsk = mean(BPSK_signal_t.^2);
fprintf('Sbpsk %d\n', Sbpsk);
% фактическая мощность зашумленного сигнала 
SNbpsk = mean(BPSK_signal_awgn.^2);
fprintf('SNbpsk %d\n', SNbpsk);
Nss = mean((BPSK_signal_t - BPSK_signal_awgn).^2);
SNR = SNbpsk/Nss;
fprintf('SNR %d\n', SNR);

В моем случае получается SNR = 1.12

Не знаю так верно или нет...

Оценку BER не делал.

Думаю BER нужно оценивать для пакета данных конкретной длины. В моем случае это 210 ms (несколько байт информации)

 

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

BER это что? Битовая ошибка?

 

https://en.wikipedia.org/wiki/Bit_error_rate

Смотрим график BER(Eb/N0), никакой полосы там нет.

 

Думаю BER нужно оценивать для пакета данных конкретной длины. В моем случае это 210 ms (несколько байт информации)

 

Для начала надо научиться просто задавать параметиры AWGN канала и проверять BER. На простейшей непрерывной модельке с известной априори синхронизацией на приёме. Никаких SNR и полос. Только BER(Eb/N0), вы должны уметь построить такой график практически, чтобы он при идеальном приёме совпадал с теоретическим.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Спасибо. В общих чертах вроде понятно.

Поскольку сталкиваюсь с этим впервые то на деле, сходу не получится.

Что такое полоса содержащая половину мощности?

1. Предполагается что приемный усилитель будет имень в своем составе полосовой фильтр с центр. частотой 36 кГц. Вопрос в том

какая у него должна быть полоса. Если заложить всего 200 Гц это будет неверно, - высокая селективность это плохая устойчивость для фазоманипулированного сигнала. Если заложить всю возможную полосу 2000 Гц которую обеспечивает гидроакустический преобразователь тоже будет не верно. Лишние спектральные компоненты в приемн. тракте тоже не нужны.

СКЗ напряжения вашего сигнала - нет представления что это ...

2. Матлаб примеры предлагают 0.00001*randn(1,N); где N кол. выборок

Что такое полоса (Fsr/2)?

3. Вроде понятно

4. Тоже понятно.

Матлаб предлагает вычислять SNR так https://www.mathworks.com/help/signal/ref/snr.html

где все завязано на частоте сигнала, кол. выборок на период, длительности сигнала.

1. Полоса пропускания предварительного фильтра приемника и эквивалентная полоса для расчета соотношения сигнал/шум это немного разные сущности.

Пример 1 - сигнал OFDM имеет прямоугольный спектр, полоса сигнала вычисляется просто dFsig=Fmax-Fmin, и никаких вопросов не возникает.

Пример 2 - сигнал имеет спектр похожий на sin(x)/x. Вот здесь и приходится обрезать спектр по каким то критериям. (Проверки для себя - спектр берут уже, проверки для сторонних - стараются обнаучить и немного спектр расширить.).

СКЗ - это среднеквадратическое значение.

2. Полоса Fsr/2 - функция randn генерирует случайный сигнал в полосе частот от 0 до половины частоты дискретизации.

0.00001*randn(1,N); - с единичным коэффициентом проще работать.

 

Матлаб предлагает вычислять SNR ... Если есть понимание физики процесса, то результат должен быть одинаков.

Я приводил расчет значений SNR по напряжению. В матлабе везде используются значения по мощности или в dB.

 

 

Вот смотрю графики BER в книжках и не вижу там никакой полосы.

 

1.Нужно измерить мощность сигнала.

2. Задать параметры блока AWGN Channel https://www.mathworks.com/help/comm/ref/awgnchannel.html

Нужный Eb/No.

Количество бит передаваемых на символ.

Мощность сигнала измеренную.

Длительность символа.

3. Практически измерить BER при заданном Eb/No и убедиться, что совпадает с теорией.

 

Eb/No = (S/N)*(W/R);

S - средняя мощности сигнала;

N - средняя мощность шума;

W - ширина полосы сигнала;

R - битовая скорость.

 

В тех расчетах, что я приводил выше SNR=sqrt((S/N), поэтому требуется использовать W в формулах напрямую.

 

 

BER это что? Битовая ошибка?

Полоса BPSK http://www.dsplib.ru/content/bpsk/bpsk.html

В матлаб SNR оцениваю так

% отношение сигнал/шум в децибелах
SNR = 0;
% смесь мод.сигнала с шумом
BPSK_signal_awgn = awgn( BPSK_signal_t, SNR);

% фактическая мощность сигнала 
Sbpsk = mean(BPSK_signal_t.^2);
fprintf('Sbpsk %d\n', Sbpsk);
% фактическая мощность зашумленного сигнала 
SNbpsk = mean(BPSK_signal_awgn.^2);
fprintf('SNbpsk %d\n', SNbpsk);
Nss = mean((BPSK_signal_t - BPSK_signal_awgn).^2);
SNR = SNbpsk/Nss;
fprintf('SNR %d\n', SNR);

В моем случае получается SNR = 1.12

Не знаю так верно или нет...

Оценку BER не делал.

Думаю BER нужно оценивать для пакета данных конкретной длины. В моем случае это 210 ms (несколько байт информации)

 

BPSK_signal_awgn = awgn( BPSK_signal_t, SNR); - Вот здесь у Вас и сидит ошибка!!!

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Eb/No = (S/N)*(W/R);

S - средняя мощности сигнала;

N - средняя мощность шума;

W - ширина полосы сигнала;

R - битовая скорость.

 

Это неправильно, в общем случае сигнал может любую полосу иметь при заданных Eb/N0, длительности символа и количестве бит в символе. Полоса - лишняя сущность, которая всё запутывает.

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

SNR  = 6:
BPSK_signal_awgn = awgn( BPSK_signal_t, SNR); - Вот здесь у Вас и сидит ошибка!!!

Но это же просто функция для подмешивания аддитивного белого гаусового шума к сигналу.

Если мощность сигнала принять за 0. Тут заложено что мощность сигнала в 2 раза больше мощности шума.

В чем может быть ошибка?

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

Вы можете написать сейчас и зарегистрироваться позже. Если у вас есть аккаунт, авторизуйтесь, чтобы опубликовать от имени своего аккаунта.

Гость
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...