acvarif 0 26 апреля, 2018 Опубликовано 26 апреля, 2018 (изменено) · Жалоба Помогите пожалуйста найти ошибку в Матлаб коде петли ФАПЧ %% Тест ФАПЧ 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 Вот что получается: Изменено 26 апреля, 2018 пользователем Acvarif Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
petrov 7 26 апреля, 2018 Опубликовано 26 апреля, 2018 · Жалоба Нарисуйте уже в симулинке вашу петлю из простейших блоков. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
acvarif 0 27 апреля, 2018 Опубликовано 27 апреля, 2018 (изменено) · Жалоба Нарисуйте уже в симулинке вашу петлю из простейших блоков. Пока смутно представляю как это сделать, поскольку в симулинк сильно плаваю. Но ошибку нашел. ФАПЧ работает, неплохо работает. Таки да, ей без разницы какой сигнал, модулированный или нет. На картинках спектры выхода NCO и входного сигнала с шумом с соотн 1, при долере 120 Гц на несущей 36 000 Гц. Остаются еще неясности. Например: 1. Как можно просто пояснить что такое натуральная частота? 2. Нет полного представления как петлю пристроить к BPSK демодулятору. Изменено 27 апреля, 2018 пользователем Acvarif Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
petrov 7 27 апреля, 2018 Опубликовано 27 апреля, 2018 · Жалоба Плавайте дальше. Вам даже картинку нарисовали в книжке. Что смутного в сумматоре, умножителе, задержке? Если с 2x2 проблема, стоит ли вообще браться за такие задачи? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
acvarif 0 27 апреля, 2018 Опубликовано 27 апреля, 2018 (изменено) · Жалоба Плавайте дальше. Вам даже картинку нарисовали в книжке. Что смутного в сумматоре, умножителе, задержке? Если с 2x2 проблема, стоит ли вообще браться за такие задачи? Не в этом дело. Любую среду нужно освоить, хотя-бы механически, сначала на простейших примерах и т. д. Что я и делаю... Вы же не будете спорить, что зная как провести один проводник к другому нельзя сразу развести плату в PCAD. Или нарисовав схему запустить ее моделирование в ModelSim. Даже неплохо владея С++ нельзя сразу писать код в Матлаб. Изменено 27 апреля, 2018 пользователем Acvarif Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
mvm54 0 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба ФАПЧ работает, неплохо работает. Таки да, ей без разницы какой сигнал, модулированный или нет. На картинках спектры выхода NCO и входного сигнала с шумом с соотн 1, при долере 120 Гц на несущей 36 000 Гц. Acvarif, Вы несколько раз приводили предельные значения SNR < 1 входного сигнала, при которых работает Ваша система. Можете привести подробный расчет ( сигнал+шум) с учетом полосы занимаемой сигналом, который использовали в моделировании. Несущая Fsig=36000 Гц; Амплитуда Asig=1 (???). Символьная Fsim=Fsig/n=100 Гц (???). (n=360 ???). Частота дискретизации Fsr=288000 Гц. ( ???) Эффективная полоса сигнала dFsig=Fsim (2Fsim ... другая ???) Гц. Формла расчета ( сигнал+шум) = .....??? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
acvarif 0 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба Acvarif, Вы несколько раз приводили предельные значения SNR < 1 входного сигнала, при которых работает Ваша система. Можете привести подробный расчет ( сигнал+шум) с учетом полосы занимаемой сигналом, который использовали в моделировании. Несущая Fsig=36000 Гц; Амплитуда Asig=1 (???). Символьная Fsim=Fsig/n=100 Гц (???). (n=360 ???). Частота дискретизации Fsr=288000 Гц. ( ???) Эффективная полоса сигнала dFsig=Fsim (2Fsim ... другая ???) Гц. Формла расчета ( сигнал+шум) = .....??? Нет. С учетом полосы занимаемой сигналом нет. Полоса BPSK сигнала равна удвоенной символьной частоте, тоесть 200 Гц В модели просто смешиваю сигнал с шумом 1:1 с амплитудой сигнала 1. Очевидно что это не совсем верно. Буду признателен если подскажете как правильно сделать расчет. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
mvm54 0 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба Нет. С учетом полосы занимаемой сигналом нет. Полоса 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)); Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
acvarif 0 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба Эффективная полоса сигнала в 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 где все завязано на частоте сигнала, кол. выборок на период, длительности сигнала. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
petrov 7 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба Эффективная полоса... Вот смотрю графики BER в книжках и не вижу там никакой полосы. 1.Нужно измерить мощность сигнала. 2. Задать параметры блока AWGN Channel https://www.mathworks.com/help/comm/ref/awgnchannel.html Нужный Eb/No. Количество бит передаваемых на символ. Мощность сигнала измеренную. Длительность символа. 3. Практически измерить BER при заданном Eb/No и убедиться, что совпадает с теорией. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
acvarif 0 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба Вот смотрю графики 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 (несколько байт информации) Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
petrov 7 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба BER это что? Битовая ошибка? https://en.wikipedia.org/wiki/Bit_error_rate Смотрим график BER(Eb/N0), никакой полосы там нет. Думаю BER нужно оценивать для пакета данных конкретной длины. В моем случае это 210 ms (несколько байт информации) Для начала надо научиться просто задавать параметиры AWGN канала и проверять BER. На простейшей непрерывной модельке с известной априори синхронизацией на приёме. Никаких SNR и полос. Только BER(Eb/N0), вы должны уметь построить такой график практически, чтобы он при идеальном приёме совпадал с теоретическим. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
mvm54 0 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба Спасибо. В общих чертах вроде понятно. Поскольку сталкиваюсь с этим впервые то на деле, сходу не получится. Что такое полоса содержащая половину мощности? 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); - Вот здесь у Вас и сидит ошибка!!! Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
petrov 7 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба Eb/No = (S/N)*(W/R); S - средняя мощности сигнала; N - средняя мощность шума; W - ширина полосы сигнала; R - битовая скорость. Это неправильно, в общем случае сигнал может любую полосу иметь при заданных Eb/N0, длительности символа и количестве бит в символе. Полоса - лишняя сущность, которая всё запутывает. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
acvarif 0 1 мая, 2018 Опубликовано 1 мая, 2018 · Жалоба SNR = 6: BPSK_signal_awgn = awgn( BPSK_signal_t, SNR); - Вот здесь у Вас и сидит ошибка!!! Но это же просто функция для подмешивания аддитивного белого гаусового шума к сигналу. Если мощность сигнала принять за 0. Тут заложено что мощность сигнала в 2 раза больше мощности шума. В чем может быть ошибка? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться