реклама на сайте
подробности

 
 
2 страниц V   1 2 >  
Reply to this topicStart new topic
> Не стабилизируется ФАПЧ в Матлаб
Acvarif
сообщение Apr 26 2018, 13:12
Сообщение #1


Знающий
****

Группа: Участник
Сообщений: 998
Регистрация: 27-08-08
Пользователь №: 39 850



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

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', cool.gif;
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
Вот что получается:
Прикрепленное изображение


Сообщение отредактировал Acvarif - Apr 26 2018, 13:20
Go to the top of the page
 
+Quote Post
petrov
сообщение Apr 26 2018, 13:31
Сообщение #2


Гуру
******

Группа: Свой
Сообщений: 2 213
Регистрация: 21-10-04
Из: Balakhna
Пользователь №: 937



Нарисуйте уже в симулинке вашу петлю из простейших блоков.
Go to the top of the page
 
+Quote Post
Acvarif
сообщение Apr 27 2018, 10:18
Сообщение #3


Знающий
****

Группа: Участник
Сообщений: 998
Регистрация: 27-08-08
Пользователь №: 39 850



Цитата(petrov @ Apr 26 2018, 16:31) *
Нарисуйте уже в симулинке вашу петлю из простейших блоков.

Пока смутно представляю как это сделать, поскольку в симулинк сильно плаваю.
Но ошибку нашел. ФАПЧ работает, неплохо работает. Таки да, ей без разницы какой сигнал, модулированный или нет.
На картинках спектры выхода NCO и входного сигнала с шумом с соотн 1, при долере 120 Гц на несущей 36 000 Гц.
Прикрепленное изображение

Остаются еще неясности.
Например:
1. Как можно просто пояснить что такое натуральная частота?
2. Нет полного представления как петлю пристроить к BPSK демодулятору.

Сообщение отредактировал Acvarif - Apr 27 2018, 10:23
Go to the top of the page
 
+Quote Post
petrov
сообщение Apr 27 2018, 10:55
Сообщение #4


Гуру
******

Группа: Свой
Сообщений: 2 213
Регистрация: 21-10-04
Из: Balakhna
Пользователь №: 937



Плавайте дальше. Вам даже картинку нарисовали в книжке. Что смутного в сумматоре, умножителе, задержке? Если с 2x2 проблема, стоит ли вообще браться за такие задачи?
Go to the top of the page
 
+Quote Post
Acvarif
сообщение Apr 27 2018, 11:48
Сообщение #5


Знающий
****

Группа: Участник
Сообщений: 998
Регистрация: 27-08-08
Пользователь №: 39 850



Цитата(petrov @ Apr 27 2018, 13:55) *
Плавайте дальше. Вам даже картинку нарисовали в книжке. Что смутного в сумматоре, умножителе, задержке? Если с 2x2 проблема, стоит ли вообще браться за такие задачи?

Не в этом дело. Любую среду нужно освоить, хотя-бы механически, сначала на простейших примерах и т. д. Что я и делаю...
Вы же не будете спорить, что зная как провести один проводник к другому нельзя сразу развести плату в PCAD.
Или нарисовав схему запустить ее моделирование в ModelSim.
Даже неплохо владея С++ нельзя сразу писать код в Матлаб.

Сообщение отредактировал Acvarif - Apr 27 2018, 11:50
Go to the top of the page
 
+Quote Post
mvm54
сообщение May 1 2018, 06:08
Сообщение #6


Участник
*

Группа: Участник
Сообщений: 40
Регистрация: 29-11-07
Пользователь №: 32 817



Цитата(Acvarif @ Apr 27 2018, 13:18) *
ФАПЧ работает, неплохо работает. Таки да, ей без разницы какой сигнал, модулированный или нет.
На картинках спектры выхода NCO и входного сигнала с шумом с соотн 1, при долере 120 Гц на несущей 36 000 Гц.


Acvarif, Вы несколько раз приводили предельные значения SNR < 1 входного сигнала, при которых работает Ваша система. Можете привести подробный расчет ( сигнал+шум) с учетом полосы занимаемой сигналом, который использовали в моделировании.
Несущая Fsig=36000 Гц; Амплитуда Asig=1 (???).
Символьная Fsim=Fsig/n=100 Гц (???). (n=360 ???).
Частота дискретизации Fsr=288000 Гц. ( ???)
Эффективная полоса сигнала dFsig=Fsim (2Fsim ... другая ???) Гц.
Формла расчета ( сигнал+шум) = .....???
Go to the top of the page
 
+Quote Post
Acvarif
сообщение May 1 2018, 06:49
Сообщение #7


Знающий
****

Группа: Участник
Сообщений: 998
Регистрация: 27-08-08
Пользователь №: 39 850



Цитата(mvm54 @ May 1 2018, 09:08) *
Acvarif, Вы несколько раз приводили предельные значения SNR < 1 входного сигнала, при которых работает Ваша система. Можете привести подробный расчет ( сигнал+шум) с учетом полосы занимаемой сигналом, который использовали в моделировании.
Несущая Fsig=36000 Гц; Амплитуда Asig=1 (???).
Символьная Fsim=Fsig/n=100 Гц (???). (n=360 ???).
Частота дискретизации Fsr=288000 Гц. ( ???)
Эффективная полоса сигнала dFsig=Fsim (2Fsim ... другая ???) Гц.
Формла расчета ( сигнал+шум) = .....???

Нет. С учетом полосы занимаемой сигналом нет. Полоса BPSK сигнала равна удвоенной символьной частоте, тоесть 200 Гц
В модели просто смешиваю сигнал с шумом 1:1 с амплитудой сигнала 1. Очевидно что это не совсем верно.
Буду признателен если подскажете как правильно сделать расчет.
Go to the top of the page
 
+Quote Post
mvm54
сообщение May 1 2018, 10:03
Сообщение #8


Участник
*

Группа: Участник
Сообщений: 40
Регистрация: 29-11-07
Пользователь №: 32 817



Цитата(Acvarif @ May 1 2018, 09:49) *
Нет. С учетом полосы занимаемой сигналом нет. Полоса 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));

Go to the top of the page
 
+Quote Post
Acvarif
сообщение May 1 2018, 10:46
Сообщение #9


Знающий
****

Группа: Участник
Сообщений: 998
Регистрация: 27-08-08
Пользователь №: 39 850



Цитата(mvm54 @ May 1 2018, 13:03) *
Эффективная полоса сигнала в 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
где все завязано на частоте сигнала, кол. выборок на период, длительности сигнала.

Go to the top of the page
 
+Quote Post
petrov
сообщение May 1 2018, 11:05
Сообщение #10


Гуру
******

Группа: Свой
Сообщений: 2 213
Регистрация: 21-10-04
Из: Balakhna
Пользователь №: 937



Цитата(mvm54 @ May 1 2018, 13:03) *
Эффективная полоса...


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

1.Нужно измерить мощность сигнала.
2. Задать параметры блока AWGN Channel https://www.mathworks.com/help/comm/ref/awgnchannel.html
Нужный Eb/No.
Количество бит передаваемых на символ.
Мощность сигнала измеренную.
Длительность символа.
3. Практически измерить BER при заданном Eb/No и убедиться, что совпадает с теорией.
Go to the top of the page
 
+Quote Post
Acvarif
сообщение May 1 2018, 12:14
Сообщение #11


Знающий
****

Группа: Участник
Сообщений: 998
Регистрация: 27-08-08
Пользователь №: 39 850



Цитата(petrov @ May 1 2018, 14:05) *
Вот смотрю графики 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 (несколько байт информации)

Go to the top of the page
 
+Quote Post
petrov
сообщение May 1 2018, 12:37
Сообщение #12


Гуру
******

Группа: Свой
Сообщений: 2 213
Регистрация: 21-10-04
Из: Balakhna
Пользователь №: 937



Цитата(Acvarif @ May 1 2018, 15:14) *
BER это что? Битовая ошибка?


https://en.wikipedia.org/wiki/Bit_error_rate
Смотрим график BER(Eb/N0), никакой полосы там нет.

Цитата(Acvarif @ May 1 2018, 15:14) *
Думаю BER нужно оценивать для пакета данных конкретной длины. В моем случае это 210 ms (несколько байт информации)


Для начала надо научиться просто задавать параметиры AWGN канала и проверять BER. На простейшей непрерывной модельке с известной априори синхронизацией на приёме. Никаких SNR и полос. Только BER(Eb/N0), вы должны уметь построить такой график практически, чтобы он при идеальном приёме совпадал с теоретическим.
Go to the top of the page
 
+Quote Post
mvm54
сообщение May 1 2018, 13:22
Сообщение #13


Участник
*

Группа: Участник
Сообщений: 40
Регистрация: 29-11-07
Пользователь №: 32 817



Цитата(Acvarif @ May 1 2018, 13:46) *
Спасибо. В общих чертах вроде понятно.
Поскольку сталкиваюсь с этим впервые то на деле, сходу не получится.
Что такое полоса содержащая половину мощности?
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.


Цитата(petrov @ May 1 2018, 14:05) *
Вот смотрю графики 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 в формулах напрямую.


Цитата(Acvarif @ May 1 2018, 15:14) *
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); - Вот здесь у Вас и сидит ошибка!!!
Go to the top of the page
 
+Quote Post
petrov
сообщение May 1 2018, 13:26
Сообщение #14


Гуру
******

Группа: Свой
Сообщений: 2 213
Регистрация: 21-10-04
Из: Balakhna
Пользователь №: 937



Цитата(mvm54 @ May 1 2018, 16:12) *
Eb/No = (S/N)*(W/R);
S - средняя мощности сигнала;
N - средняя мощность шума;
W - ширина полосы сигнала;
R - битовая скорость.


Это неправильно, в общем случае сигнал может любую полосу иметь при заданных Eb/N0, длительности символа и количестве бит в символе. Полоса - лишняя сущность, которая всё запутывает.
Go to the top of the page
 
+Quote Post
Acvarif
сообщение May 1 2018, 15:27
Сообщение #15


Знающий
****

Группа: Участник
Сообщений: 998
Регистрация: 27-08-08
Пользователь №: 39 850



Цитата(mvm54 @ May 1 2018, 16:22) *
Код
SNR  = 6:
BPSK_signal_awgn = awgn( BPSK_signal_t, SNR); - Вот здесь у Вас и сидит ошибка!!!

Но это же просто функция для подмешивания аддитивного белого гаусового шума к сигналу.
Если мощность сигнала принять за 0. Тут заложено что мощность сигнала в 2 раза больше мощности шума.
В чем может быть ошибка?
Go to the top of the page
 
+Quote Post

2 страниц V   1 2 >
Reply to this topicStart new topic
3 чел. читают эту тему (гостей: 3, скрытых пользователей: 0)
Пользователей: 0

 


RSS Текстовая версия Сейчас: 19th August 2018 - 06:43
Рейтинг@Mail.ru


Страница сгенерированна за 0.00944 секунд с 7
ELECTRONIX ©2004-2016