Моделирование ФАПЧ

abraziv
Сообщения: 48
Зарегистрирован: 08 апр 2015, 15:16

Моделирование ФАПЧ

Сообщение abraziv » 24 фев 2017, 09:44

Здравствуйте. Хотел спросить знающих людей. В результате моделирования петли ФАПЧ получил следующие диаграммы:
Снимок3.JPG
На первом плоте представлен исходный сигнал и сигнал NCO.
На втором графике представлен сигнал ошибки на выходе петлевого фильтра.
Вроде бы всё хорошо, но по фазе синхронизации почему то нет. Проблема в том, что сигнал ошибки не стабилизируются около нуля, а немного ниже. С чем это может быть связано?
for k = 2:N

nco(k) = sin(2*pi*Fnco*Ts*(k-1) + pnco(k-1));

e(k) = Kd * s(k-1) * nco(k);

phi_output(k) = (Kp+Ki*Ts/2)*e(k) + (Ki*Ts/2-Kp)*e(k-1) + phi_output(k-1);

pnco(k) = pnco(k-1) + Ko * phi_output(k-1);

end;

Хотел так же спросить как задать время вхождения в синхронизм по частоте? Есть формулы, но как их связать с коэф усиления ГУНа и ФД, ведь по логике как раз это усиление и управляет временем вхождения в синхронизм.

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 892
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Моделирование ФАПЧ

Сообщение Бахурин Сергей » 24 фев 2017, 12:36

Ну если по фазе синхронизации нет то значит все плохо. Мне кажется вы как-то странно моделируете контур. Я так понял вы пытаетесь смоделировать аналоговый контур (не цифровой), только почему то используете разностные уравнения. Боюсь что так может не прокатить. Попробуйте пример отсюда. Там правда на си написано но на матлаб переписать думаю не сложно.

abraziv
Сообщения: 48
Зарегистрирован: 08 апр 2015, 15:16

Re: Моделирование ФАПЧ

Сообщение abraziv » 24 фев 2017, 13:10

Спасибо за ответ. Пытаюсь смоделировать цифровой контур. Я взял аналоговый фильтр прототип с передаточной характеристикой: H(s) = K1 + K2/s - это активный фильтр с бесконечным усилением. После перевёл в Z область с помощью билинейного преобразования, в итоге получил:
out(k) = (K1+K2*Ts/2)*in(k) + (K2*Ts/2-K1)*in(k-1) + out(k-1)
Шумовую полосу определил путём задания временем синхронизации по частоте Tf и максимальной частотной отстройкой dF0max
B = (4*dF0max^2 / Tf)^(1/3),
где
Tf = 4*dF0max^2/(B^3)
Так корректно делать?
Дальше определил K1 и K2 по книге RICEа:
Снимок.JPG
Вот весь код:

close all;

Kd=1;
Ko=1;

Fs = 1E+6;
Ts = 1/Fs;
demp = 1/sqrt(2);
dF0max = 1000;
Tf=1000e-6;

B = (4*dF0max^2 / Tf)^(1/3);
Wn = 2*B/(demp+1/(4*demp));

df = 2*pi*sqrt(1)*demp*B

Kp = (4*demp*Wn/(1+2*demp*Wn + Wn^2))/Kd/Ko;
Ki = (4*Wn^2/(1+2*demp*Wn + Wn^2))/Kd/Ko;

N = 100000;
F0 = 1E+3;
pho = 45/180*pi;
s = sin(2 * pi * F0 / Fs * (1:N) + pho);
Fnco = F0 + 500;

e(1) = 0;
phi_output(1) = 0;
nco(1) = 0;
pnco(1) = 0;

for k = 2:N

nco(k) = sin(2*pi*Fnco*Ts*(k-1) + pnco(k-1));

e(k) = Kd * s(k-1) * nco(k);

phi_output(k) = (Kp+Ki*Ts/2)*e(k) + (Ki*Ts/2-Kp)*e(k-1) + phi_output(k-1);

pnco(k) = pnco(k-1) + Ko * phi_output(k);

end;

n = N-50000;

plot((1:n)*Ts, s(1:n),(1:n)*Ts, nco(1:n));
figure(2);
plot((1:n)*Ts, phi_output(1:n));

Вашу статью читал и код смотрел не раз, не помогло ((((((

abraziv
Сообщения: 48
Зарегистрирован: 08 апр 2015, 15:16

Re: Моделирование ФАПЧ

Сообщение abraziv » 24 фев 2017, 16:53

Нашёл другой источник, попробовал смоделировать по нему. Заработало, вот только при одном жёстко условии, а именно при Fs в 4 раза больше максимальной частоты входного сигнала (Это же бред согласитесь?), меньше и больше не работает. Скоро башка уже лопнет. Так же время синхронизации не соответствует расчётному((((((
Снимок.JPG
close all;

Kd=1;
Ko=1;

Fs = 4E+3;
Ts = 1/Fs;

demp = 1/sqrt(2);
dF0max = 1000;
Tf=100e-6;

B = (4*dF0max^2 / Tf)^(1/3)
Wn = 2*B/(demp+1/(4*demp));

a0 = 4*(Wn*Ts)^2 +demp*Wn*Ts+4;
a1 = 2*(Wn*Ts)^2 - 8;
a2 = (Wn*Ts)^2 - 4*demp*Wn*Ts+4;

go = 1/(Kd*Ko)*(a1/a0+2);
g1 = 1/(Kd*Ko)*(a2/a0-1);

N = 1000;
F0 = 1E+3;
pho = 90/180*pi;
s = sin(2 * pi * F0 * Ts * (1:N) + pho);
Fnco = F0 + 400;

e(1) = 0;
phi_output(1) = 0;
nco(1) = 0;
pnco(1) = 0;

for k = 2:N

nco(k) = sin(2*pi*Fnco*Ts*(k-1) + pnco(k-1));

e(k) = Kd * s(k-1) * nco(k);

phi_output(k) = go*e(k) + g1*e(k-1) + phi_output(k-1);

pnco(k) = Ko * phi_output(k) + pnco(k-1);

end;

n = N;

plot((1:n)*Ts, s(1:n),(1:n)*Ts, nco(1:n));
figure(2);
plot((1:n)*Ts, phi_output(1:n));

abraziv
Сообщения: 48
Зарегистрирован: 08 апр 2015, 15:16

Re: Моделирование ФАПЧ

Сообщение abraziv » 24 фев 2017, 18:34

По логике Wn << 2*pi*Fs, т.е. чем выше Fs тем лучше, что вы и писали в своей статье. Не знаю что и делать, вроде бы много литературы на этот счёт, сам вывожу эти формулы, а в итоге ничего не работает. Из-за этого не могу двигаться дальше.

Аватара пользователя
quato_a
Сообщения: 3
Зарегистрирован: 24 фев 2017, 23:57
Откуда: Зеленоград

Re: Моделирование ФАПЧ

Сообщение quato_a » 25 фев 2017, 00:00

abraziv писал(а):
24 фев 2017, 16:53
phi_output(k) = go*e(k) + g1*e(k-1) + phi_output(k-1);
phi_output(k) = go*e(k) + (g1-g0)*e(k-1) + phi_output(k-1);
Суббота начинается в понедельник

abraziv
Сообщения: 48
Зарегистрирован: 08 апр 2015, 15:16

Re: Моделирование ФАПЧ

Сообщение abraziv » 25 фев 2017, 10:56

Спасибо за ответ. Хоть кто-то отозвался. Я всё правильно написал, передаточная характеристика фильтра F(z) = (g0+g1*z^(-1))/(1-z^(-1))
https://theses.lib.vt.edu/theses/availa ... ed/etd.pdf
18 страница.
Если моделировать как
phi_output(k) = go*e(k) + (g1-g0)*e(k-1) + phi_output(k-1);
то вообще не работает.

abraziv
Сообщения: 48
Зарегистрирован: 08 апр 2015, 15:16

Re: Моделирование ФАПЧ

Сообщение abraziv » 25 фев 2017, 11:07

Подвох здесь, как-то связан с резонансной частотой Wn, только не могу понять как.

Аватара пользователя
quato_a
Сообщения: 3
Зарегистрирован: 24 фев 2017, 23:57
Откуда: Зеленоград

Re: Моделирование ФАПЧ

Сообщение quato_a » 25 фев 2017, 12:54

У меня тоже поначалу работа петлевого фильтра была не устойчива. Потом рассчитал коэффициенты по http://www.dsplib.ru/content/dpll/dpll.html, которая здесь уже указывалась. Заметил, что при увеличении/уменьшении резонансной частоты повышается/понижается скорость сходимости по частоте. Но при увеличении скорости поиска/слежения повышался шум в полной фазе NCO (в конце: http://sernam.ru/book_p_net.php?id=101).

Сейчас успешно работает в железе в схеме Костоса при fs =~ 75 МГц, f0 =~ 6.3 МГц и df = +/-300 кГц.

Код: Выделить всё

% dsplib.ru DPLL example

% close all;
clear all;
clc;

Ks = 1000;

Kd = 1;
Ko = 1;
wp = 2*pi*50;
zeta = 0.5;

fs = 1e4;
T = 1/fs;
f0 = 1e3;
phi = -1.5;
df = 100;
fg = f0 - df;

teta = df/fs

g1 = 2 * (1 - exp(-wp * zeta * T) * cos(wp * sqrt(1 - zeta^2)*T));
g2 = g1 + exp(-2 * wp * zeta * T) - 1;

Ki = g2 / (Ko*Kd);
Kp = g1 / (Ko*Kd);

N = 3000;

p = zeros(1,N); % фаза NCO
v = zeros(1,N); % выход фазового детектора
e = zeros(1,N); % выход петлевого фильтра
y = zeros(1,N); % сигнал на выходе NCO
r = zeros(1,N); % ошибка слежения

% исходный сигнал
s = Ks * sin(2 * pi * f0/fs * [1:1:N]);

for i=2:N
    
    % выход фазового детектора
    v(i) = Kd * s(i) * cos(2 * pi * fg/fs * i + p(i-1));
    
    % разностное уравнение петлевого фильтра
    e(i) = Kp * v(i) + (Ki-Kp) * v(i-1) + e(i-1);
    
    % разностное уравнение NCO
    p(i) = 1/Ks * Ko * e(i-1) + p(i-1);
    
    % сигнал, синхронизованный со входным
    y(i) = Ks * sin(2 * pi * fg/fs * i + p(i-1));
    
    % ошибка слежения
    r(i) = s(i) - y(i);
end

figure;
subplot 311; plot(v); grid on; title('v');
subplot 312; plot(e); grid on; title('e');
subplot 313; plot(p); grid on; title('e');

figure;
subplot 211; plot(p); grid on; hold on;
             plot(2*pi*df/fs*[1:1:N], 'r'); legend('recovery', 'ideal');
subplot 212; plot(s); grid on; hold on;
             plot(y, 'r');
             plot(r, 'g'); legend('input', 'output', 'error');
Последний раз редактировалось quato_a 25 фев 2017, 16:23, всего редактировалось 1 раз.
Суббота начинается в понедельник

abraziv
Сообщения: 48
Зарегистрирован: 08 апр 2015, 15:16

Re: Моделирование ФАПЧ

Сообщение abraziv » 25 фев 2017, 15:36

Спасибо огромное за ответ и код. Я уже пробовал моделировать согласно статьи Сергея. Не получается связать время синхронизации по частоте с результатом моделирования.

demp = 1/sqrt(2);
dF0max = 200;
Tf=100e-3;
B = (4*dF0max^2 / Tf)^(1/3);
wp = 2*B/(demp+1/(4*demp));

А так всё работает. Просто не хочется методом тыка подгонять коэффициенты фильтра.

Ответить

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 1 гость