Версия для печати темы
Нажмите сюда для просмотра этой темы в обычном формате
Форум разработчиков электроники ELECTRONIX.ru _ Вопросы системного уровня проектирования _ BPSK матлаб модулятор демодулятор
Автор: Acvarif Mar 19 2018, 12:14
Имеется код BPSK модема, передатчик и приемник.
CODE
clear all;
Br = 117.1875*10^3;
Fd = Br*8;
Fif = 20*10^6;
Fdif = 60*10^6;
Fdopl = 100;
%% === TRANSMITTER ===
M = pngen(5);
D = rand(1, 200*8) > 0.5;
CRC = zeros(1, 16);
for i=1:length(D),
CRC(16) = xor(CRC(16), D(i));
CRC = [CRC(2:16) CRC(1)];
end;
FILL = ones(1, 200)*0.5;
DATA = [dec2bin(length(D)/8, 11)=='1' D CRC];
packet = [FILL M M diff_enc([0 DATA]) FILL]*2-1;
% redescr
Ndkr = Fd/Br;
pkt = packet(ceil(0.01:1/Ndkr:length(packet)));
hch = firls(31, [0 .24 .3 .5]*2, [1 1 0 0]);
pktf = filter(hch, [1], pkt);
% upsample
t1 = upsample4(pktf, 31);
t2 = upsample4(t1, 31);
pktfu = upsample4(t2, 31)';
% -> IF
l = length(pktfu);
s = sin(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
c = cos(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
pkts = pktfu.*s(1:l);
pktc = pktfu.*c(1:l);
%% === TRANSMIT LINE ===
Ir = awgn(awgn(pkts, 0.05), 0.01);
Qr = awgn(awgn(pktc, 0.05), 0.01);
%% === RECEIVER ===
% -> 0
l = length(Ir);
s = sin(2*pi*Fif*(0:1/Fdif:l/Fdif));
c = cos(2*pi*Fif*(0:1/Fdif:l/Fdif));
Rreu = Ir.*s(1:l) + Qr.*c(1:l);
Rimu = Qr.*s(1:l) - Ir.*c(1:l);
tr00 = downsample4(Rreu);
tr01 = downsample4(tr00);
Rre_ = downsample4(tr01);
Rre = filter(hch, [1], Rre_);
tr10 = downsample4(Rimu);
tr11 = downsample4(tr10);
Rim_ = downsample4(tr11);
Rim = filter(hch, [1], Rim_);
%figure; hold on; plot(Rre); plot(Rim, 'r');
% detection
Ms = M(ceil(0.01:1/Ndkr:length(M)));
cfh = Ms(length(Ms):-1:1);
cre = filter(cfh, [1], Rre)/length(cfh);
cim = filter(cfh, [1], Rim)/length(cfh);
%cr = sqrt(cre.^2 + cim.^2);
cr = (7/8)*max(abs(cre), abs(cim)) + (1/8)*min(abs(cre), abs(cim));
%figure; hold on; plot(cre); plot(cim); plot(cr, 'r');
env(length(cr)) = 0;
for j=1:length(cr),
f = max(1, j-Ndkr*length(M)+2*Ndkr);
t = min(length(cr), max(2, j-Ndkr));
env(j) = sum(cr(f:t))/(t-f);
end;
%figure; hold on; plot(cr, 'r'); plot(env, 'g'); plot(env*5);
%
for i=length(M)*Ndkr+1:length(cr)-1,
c1 = cr(i) > env(i)*5;
% c2 = (cr(i) > cr(i-1)) && (cr(i) > cr(i+1));
c2 = (cr(i-1) > cr(i-2)) && (cr(i-1) > cr(i));
c3 = abs(cr(i)-cr(i-length(M)*Ndkr)) < 4*env(i);
c4 = cr(i-length(M)*Ndkr) > env(i-length(M)*Ndkr)*5;
if c1 && c2 && c3 && c4,
break;
end;
end;
i
re = Rre(i:length(Rre));
im = Rim(i:length(Rim));
%figure; hold on; plot(re); plot(im, 'r');
figure; hold on; plot(cr, 'r'); plot(env, 'g'); plot(env*5); plot((c1&c2&c3&c4)*1.3, 'k');
%
l = length(re);
lb = length(packet);
Yre(lb) = 0;
Yim(lb) = 0;
fre(lb) = 0;
fim(lb) = 0;
are(lb) = 0;
aim(lb) = 0;
mark(l) = 0;
bit(lb) = 0;
k = 1;
gamma = 0.1;
Yre(1) = re(1);
Yim(1) = im(1);
for i = 8:8:length(re)-10,
mark(i) = 1;
are(k) = sum(re(i-6:i-1));
aim(k) = sum(im(i-6:i-1));
% okre(k) = are(k)*Yim(k) - aim(k)*Yre(k);
okim(k) = aim(k)*Yim(k) + are(k)*Yre(k);
bit(k) = sign(okim(k));
fre(k) = are(k)*bit(k);
fim(k) = aim(k)*bit(k);
Yre(k+1) = (1-gamma)*Yre(k) + gamma*fre(k);
Yim(k+1) = (1-gamma)*Yim(k) + gamma*fim(k);
k = k + 1;
end;
bit_dec = diff_dec(sign(bit+1));
bits = bit_dec(2:length(bit_dec));
%figure; hold on; plot(okim); plot(okre, 'r');
%figure; hold on; plot(re); plot(Yre(ceil(0.01:1/Ndkr:length(bit))), 'r');
%figure; hold on; plot(im); plot(Yim(ceil(0.01:1/Ndkr:length(bit))), 'r');
figure; hold on; plot(im); plot(mark, 'k');
%
len = bin2dec(int2str(bits(1:11)));
data = bits(12:12+len*8-1);
rcrc = bits(12+len*8:12+len*8+15);
crc = zeros(1, 16);
for i=1:length(D),
crc(16) = xor(crc(16), data(i));
crc = [crc(2:16) crc(1)];
end;
disp('---');
disp(sprintf('Received packet length: %d bytes', len));
disp(sprintf('Data ok: %d', sum(abs(data-D))==0));
disp(sprintf('CRC ok: %d', sum(abs(crc-rcrc))==0));
Br = 117.1875*10^3;
Fd = Br*8;
Fif = 20*10^6;
Fdif = 60*10^6;
Fdopl = 100;
%% === TRANSMITTER ===
M = pngen(5);
D = rand(1, 200*8) > 0.5;
CRC = zeros(1, 16);
for i=1:length(D),
CRC(16) = xor(CRC(16), D(i));
CRC = [CRC(2:16) CRC(1)];
end;
FILL = ones(1, 200)*0.5;
DATA = [dec2bin(length(D)/8, 11)=='1' D CRC];
packet = [FILL M M diff_enc([0 DATA]) FILL]*2-1;
% redescr
Ndkr = Fd/Br;
pkt = packet(ceil(0.01:1/Ndkr:length(packet)));
hch = firls(31, [0 .24 .3 .5]*2, [1 1 0 0]);
pktf = filter(hch, [1], pkt);
% upsample
t1 = upsample4(pktf, 31);
t2 = upsample4(t1, 31);
pktfu = upsample4(t2, 31)';
% -> IF
l = length(pktfu);
s = sin(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
c = cos(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
pkts = pktfu.*s(1:l);
pktc = pktfu.*c(1:l);
%% === TRANSMIT LINE ===
Ir = awgn(awgn(pkts, 0.05), 0.01);
Qr = awgn(awgn(pktc, 0.05), 0.01);
%% === RECEIVER ===
% -> 0
l = length(Ir);
s = sin(2*pi*Fif*(0:1/Fdif:l/Fdif));
c = cos(2*pi*Fif*(0:1/Fdif:l/Fdif));
Rreu = Ir.*s(1:l) + Qr.*c(1:l);
Rimu = Qr.*s(1:l) - Ir.*c(1:l);
tr00 = downsample4(Rreu);
tr01 = downsample4(tr00);
Rre_ = downsample4(tr01);
Rre = filter(hch, [1], Rre_);
tr10 = downsample4(Rimu);
tr11 = downsample4(tr10);
Rim_ = downsample4(tr11);
Rim = filter(hch, [1], Rim_);
%figure; hold on; plot(Rre); plot(Rim, 'r');
% detection
Ms = M(ceil(0.01:1/Ndkr:length(M)));
cfh = Ms(length(Ms):-1:1);
cre = filter(cfh, [1], Rre)/length(cfh);
cim = filter(cfh, [1], Rim)/length(cfh);
%cr = sqrt(cre.^2 + cim.^2);
cr = (7/8)*max(abs(cre), abs(cim)) + (1/8)*min(abs(cre), abs(cim));
%figure; hold on; plot(cre); plot(cim); plot(cr, 'r');
env(length(cr)) = 0;
for j=1:length(cr),
f = max(1, j-Ndkr*length(M)+2*Ndkr);
t = min(length(cr), max(2, j-Ndkr));
env(j) = sum(cr(f:t))/(t-f);
end;
%figure; hold on; plot(cr, 'r'); plot(env, 'g'); plot(env*5);
%
for i=length(M)*Ndkr+1:length(cr)-1,
c1 = cr(i) > env(i)*5;
% c2 = (cr(i) > cr(i-1)) && (cr(i) > cr(i+1));
c2 = (cr(i-1) > cr(i-2)) && (cr(i-1) > cr(i));
c3 = abs(cr(i)-cr(i-length(M)*Ndkr)) < 4*env(i);
c4 = cr(i-length(M)*Ndkr) > env(i-length(M)*Ndkr)*5;
if c1 && c2 && c3 && c4,
break;
end;
end;
i
re = Rre(i:length(Rre));
im = Rim(i:length(Rim));
%figure; hold on; plot(re); plot(im, 'r');
figure; hold on; plot(cr, 'r'); plot(env, 'g'); plot(env*5); plot((c1&c2&c3&c4)*1.3, 'k');
%
l = length(re);
lb = length(packet);
Yre(lb) = 0;
Yim(lb) = 0;
fre(lb) = 0;
fim(lb) = 0;
are(lb) = 0;
aim(lb) = 0;
mark(l) = 0;
bit(lb) = 0;
k = 1;
gamma = 0.1;
Yre(1) = re(1);
Yim(1) = im(1);
for i = 8:8:length(re)-10,
mark(i) = 1;
are(k) = sum(re(i-6:i-1));
aim(k) = sum(im(i-6:i-1));
% okre(k) = are(k)*Yim(k) - aim(k)*Yre(k);
okim(k) = aim(k)*Yim(k) + are(k)*Yre(k);
bit(k) = sign(okim(k));
fre(k) = are(k)*bit(k);
fim(k) = aim(k)*bit(k);
Yre(k+1) = (1-gamma)*Yre(k) + gamma*fre(k);
Yim(k+1) = (1-gamma)*Yim(k) + gamma*fim(k);
k = k + 1;
end;
bit_dec = diff_dec(sign(bit+1));
bits = bit_dec(2:length(bit_dec));
%figure; hold on; plot(okim); plot(okre, 'r');
%figure; hold on; plot(re); plot(Yre(ceil(0.01:1/Ndkr:length(bit))), 'r');
%figure; hold on; plot(im); plot(Yim(ceil(0.01:1/Ndkr:length(bit))), 'r');
figure; hold on; plot(im); plot(mark, 'k');
%
len = bin2dec(int2str(bits(1:11)));
data = bits(12:12+len*8-1);
rcrc = bits(12+len*8:12+len*8+15);
crc = zeros(1, 16);
for i=1:length(D),
crc(16) = xor(crc(16), data(i));
crc = [crc(2:16) crc(1)];
end;
disp('---');
disp(sprintf('Received packet length: %d bytes', len));
disp(sprintf('Data ok: %d', sum(abs(data-D))==0));
disp(sprintf('CRC ok: %d', sum(abs(crc-rcrc))==0));
При компиляции возникает ошибка:
Undefined function 'pngen' for input arguments of type 'double'.
Error in Modem (line 10)
M = pngen(5);
Я так понимаю, что функция pngen должна быть где-то описана. Или такая функция уже встроена в Матлаб последних версий.
Использую R2014a.
Подскажите пожалуйста в чем проблема.
Если кто сталкивался, подсобите пожалуйста простейшим рабочим кодом BPSK модема (передатчик + приемник).
Автор: Grizzzly Mar 19 2018, 12:18
Нет, это не встроенная функция. Ее написал автор этого кода. Надо искать реализацию.
Автор: Acvarif Mar 19 2018, 13:00
Цитата(Grizzzly @ Mar 19 2018, 16:18) 
Нет, это не встроенная функция. Ее написал автор этого кода. Надо искать реализацию.
Спасибо. Понятно.
CODE
clear
% несущаЯ
Fc = 2400;
% частота выборки бит из вектора Х
Fd = 2400;
% частота выборки выходного модулированного сигнала
Fs = 19200;
% количество сэмплов на 1 бит
Sb = Fs/Fd;
% количество 0 сэмплов перед сигналом
ZeroS = 10;%round(Sb/0.8);
% количество 0 сэмплов перед сигналом
NoiseLev = 0.2;
% ограничение тракта
Amax = 5;
% Количество наборов фаз на символ
M = 2;
% МодулирующаЯ последовательность
x = repmat([0 1 0 0 0 0 0 1 1 1 1 1 0 1 0 1 0 1 0 0 1 0 0 0],1,1); x=x(
;
% Modulate, keeping track of time.
[y,t] = dmod(x,Fc,Fd,[Fs pi/2],'psk',M);
% Всего выборок в выходом сигнале
samples_total = length(y);
if ZeroS > 0
y(ZeroS + 1: samples_total) = y (1 : samples_total - ZeroS);
y(1 : ZeroS) = 0;
end
noise = NoiseLev*randn(size(y));
% yn сигнал на входе демодулЯтора
yn = y + noise;
% Смоделировать выход за динамический диапазон по входу
for sn = 1 : samples_total
if yn(sn) > Amax
yn(sn) = Amax;
elseif yn(sn) < -Amax
yn(sn) = -Amax;
end
end
% ДемодулЯциЯ
% Вектор номеров выборок
samples_vector = 1 : samples_total;
%ГенерациЯ векторов sin и cos локальной несущей длЯ коррелЯции c Y
dFs = 2*pi*Fc/Fs;
% sin и cos на длине одного периода
sinFc = sin(0:dFs:2*pi-dFs)';
cosFc = cos(0:dFs:2*pi-dFs)';
samples_sin = length(sinFc);
% sin и cos на длине вектора входного сигнала
sinY = repmat(sinFc, samples_total/samples_sin, 1);
cosY = repmat(cosFc, samples_total/samples_sin, 1);
% ДвигаЯсь по выборкам скользЯщим окном, равным длине битового интервала, вычислЯем коррелЯцию с несущей
demod_y(1:samples_total) = 0;
demod_I(1:samples_total) = yn.*cosY;
demod_Q(1:samples_total) = yn.*sinY;
demod_y = demod_I + demod_Q*j;
% Усреднение I и Q на длине битового интервала
demod_avr_I(1:samples_total) = 0;
demod_avr_Q(1:samples_total) = 0;
avr_interval = samples_sin;
for sn = avr_interval : samples_total
demod_avr_I(sn) = sum(demod_I( sn - avr_interval + 1: sn ))/avr_interval*2;
demod_avr_Q(sn) = sum(demod_Q( sn - avr_interval + 1: sn ))/avr_interval*2;
end
% Перемножение усредненных I и Q с задержанной на 1 бит версией
demod_dly_I(1:samples_total) = 0;
demod_dly_Q(1:samples_total) = 0;
demod_dly_I(1+Sb:samples_total) = demod_avr_I(1:samples_total - Sb).* demod_avr_I(1+Sb:samples_total);
demod_dly_Q(1+Sb:samples_total) = demod_avr_Q(1:samples_total - Sb).* demod_avr_Q(1+Sb:samples_total);
demod_dly = demod_dly_I + demod_dly_Q;
% ФильтрациЯ демодулированного сигнала
[demod_b,demod_a]=butter(2,(Fd/Fs));
demod_filt = filter(demod_b,demod_a,sign(demod_dly));
subplot(5,1,1), plot (samples_vector, y, 'b', samples_vector, y, 'r'), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,2), plot (samples_vector, yn), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,3), plot (samples_vector, demod_I, samples_vector, demod_Q, samples_vector, abs(demod_y)), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,4), plot (samples_vector, demod_dly_I, samples_vector, demod_dly_Q, samples_vector, demod_dly), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,5), plot (samples_vector, demod_dly, samples_vector, sign(demod_dly), samples_vector, sign(demod_filt)), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]',...
'YLimMode', 'manual',...
'YLim', [-1.5, 1.5]',...
'YTickMode', 'manual',...
'YTick', [-1.5:0.5:1.5]');
% несущаЯ
Fc = 2400;
% частота выборки бит из вектора Х
Fd = 2400;
% частота выборки выходного модулированного сигнала
Fs = 19200;
% количество сэмплов на 1 бит
Sb = Fs/Fd;
% количество 0 сэмплов перед сигналом
ZeroS = 10;%round(Sb/0.8);
% количество 0 сэмплов перед сигналом
NoiseLev = 0.2;
% ограничение тракта
Amax = 5;
% Количество наборов фаз на символ
M = 2;
% МодулирующаЯ последовательность
x = repmat([0 1 0 0 0 0 0 1 1 1 1 1 0 1 0 1 0 1 0 0 1 0 0 0],1,1); x=x(
;% Modulate, keeping track of time.
[y,t] = dmod(x,Fc,Fd,[Fs pi/2],'psk',M);
% Всего выборок в выходом сигнале
samples_total = length(y);
if ZeroS > 0
y(ZeroS + 1: samples_total) = y (1 : samples_total - ZeroS);
y(1 : ZeroS) = 0;
end
noise = NoiseLev*randn(size(y));
% yn сигнал на входе демодулЯтора
yn = y + noise;
% Смоделировать выход за динамический диапазон по входу
for sn = 1 : samples_total
if yn(sn) > Amax
yn(sn) = Amax;
elseif yn(sn) < -Amax
yn(sn) = -Amax;
end
end
% ДемодулЯциЯ
% Вектор номеров выборок
samples_vector = 1 : samples_total;
%ГенерациЯ векторов sin и cos локальной несущей длЯ коррелЯции c Y
dFs = 2*pi*Fc/Fs;
% sin и cos на длине одного периода
sinFc = sin(0:dFs:2*pi-dFs)';
cosFc = cos(0:dFs:2*pi-dFs)';
samples_sin = length(sinFc);
% sin и cos на длине вектора входного сигнала
sinY = repmat(sinFc, samples_total/samples_sin, 1);
cosY = repmat(cosFc, samples_total/samples_sin, 1);
% ДвигаЯсь по выборкам скользЯщим окном, равным длине битового интервала, вычислЯем коррелЯцию с несущей
demod_y(1:samples_total) = 0;
demod_I(1:samples_total) = yn.*cosY;
demod_Q(1:samples_total) = yn.*sinY;
demod_y = demod_I + demod_Q*j;
% Усреднение I и Q на длине битового интервала
demod_avr_I(1:samples_total) = 0;
demod_avr_Q(1:samples_total) = 0;
avr_interval = samples_sin;
for sn = avr_interval : samples_total
demod_avr_I(sn) = sum(demod_I( sn - avr_interval + 1: sn ))/avr_interval*2;
demod_avr_Q(sn) = sum(demod_Q( sn - avr_interval + 1: sn ))/avr_interval*2;
end
% Перемножение усредненных I и Q с задержанной на 1 бит версией
demod_dly_I(1:samples_total) = 0;
demod_dly_Q(1:samples_total) = 0;
demod_dly_I(1+Sb:samples_total) = demod_avr_I(1:samples_total - Sb).* demod_avr_I(1+Sb:samples_total);
demod_dly_Q(1+Sb:samples_total) = demod_avr_Q(1:samples_total - Sb).* demod_avr_Q(1+Sb:samples_total);
demod_dly = demod_dly_I + demod_dly_Q;
% ФильтрациЯ демодулированного сигнала
[demod_b,demod_a]=butter(2,(Fd/Fs));
demod_filt = filter(demod_b,demod_a,sign(demod_dly));
subplot(5,1,1), plot (samples_vector, y, 'b', samples_vector, y, 'r'), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,2), plot (samples_vector, yn), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,3), plot (samples_vector, demod_I, samples_vector, demod_Q, samples_vector, abs(demod_y)), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,4), plot (samples_vector, demod_dly_I, samples_vector, demod_dly_Q, samples_vector, demod_dly), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,5), plot (samples_vector, demod_dly, samples_vector, sign(demod_dly), samples_vector, sign(demod_filt)), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]',...
'YLimMode', 'manual',...
'YLim', [-1.5, 1.5]',...
'YTickMode', 'manual',...
'YTick', [-1.5:0.5:1.5]');
Очевидно в этом коде функция dmod тоже принадлежит автору?
http://electronix.ru/redirect.php?https://www.gaussianwaves.com/digital-modulations-using-matlab-build-simulation-models-from-scratch/
Я так понимаю это за деньги. Можно-ли найти подобную литературу не за деньги?
Хотя dmod вроде существует http://electronix.ru/redirect.php?http://matlab.exponenta.ru/communication/book2/6/dmod.php
Но начиная с какой версии Матлаб?
Даже код кое-какой имеется http://electronix.ru/redirect.php?http://read.pudn.com/downloads40/sourcecode/book/139496/comm/COMM/DMOD.M__.htm
Автор: ataradov Mar 19 2018, 16:59
Ух ты. Я написал этот код много лет назад и все эти функции есть у меня на домашнем компе.
Я так предполагаю это из РТИ? Или это код дальше разошелся?
pngen() генерирует М-последовательность длинны 2^n, где n - это параметр функции. Можно просто найти любую М-последовательность длинны 32 и подставить напрямую.
Автор: ataradov Mar 20 2018, 00:45
Вот исходник этой pngen():
CODE
function x=pngen(p);
% x=pnegn(p) - Generate pseudo-random sequence
% p = either a scaler (2..12) to generate a PN sequence
% using a predefined polynomial of length p, or
% a binary generating polynomial.
% x = sequence, of length 2^length(p) - 1.
% Predefined seqence polynomials:
pl{2} = [1 1];
pl{3} = [1 1 0];
pl{4} = [1 0 0 1];
pl{5} = [1 0 1 1 1];
pl{6} = [1 1 1 0 0 1];
pl{7} = [1 1 0 0 1 0 1];
pl{8} = [1 0 1 0 1 1 1 1];
pl{9} = [1 0 1 0 1 1 0 1 1];
pl{10} = [1 0 1 0 1 0 1 0 1 1];
pl{11} = [1 0 1 0 1 0 1 0 1 0 1];
pl{12} = [1 1 0 1 0 0 0 0 1 1 0 1];
if length(p) == 1
p = pl{p};
end
L=length(p);
N=2^L-1;
s=ones(1,length(p));
x=zeros(1,N);
for i=1:N
snew = mod(sum(and(p,s)),2);
x(i) = snew;
s=[s(2:end),snew];
end
% x=pnegn(p) - Generate pseudo-random sequence
% p = either a scaler (2..12) to generate a PN sequence
% using a predefined polynomial of length p, or
% a binary generating polynomial.
% x = sequence, of length 2^length(p) - 1.
% Predefined seqence polynomials:
pl{2} = [1 1];
pl{3} = [1 1 0];
pl{4} = [1 0 0 1];
pl{5} = [1 0 1 1 1];
pl{6} = [1 1 1 0 0 1];
pl{7} = [1 1 0 0 1 0 1];
pl{8} = [1 0 1 0 1 1 1 1];
pl{9} = [1 0 1 0 1 1 0 1 1];
pl{10} = [1 0 1 0 1 0 1 0 1 1];
pl{11} = [1 0 1 0 1 0 1 0 1 0 1];
pl{12} = [1 1 0 1 0 0 0 0 1 1 0 1];
if length(p) == 1
p = pl{p};
end
L=length(p);
N=2^L-1;
s=ones(1,length(p));
x=zeros(1,N);
for i=1:N
snew = mod(sum(and(p,s)),2);
x(i) = snew;
s=[s(2:end),snew];
end
Если нужна помощь с остальным - обращайтесь, помогу.
Автор: Acvarif Mar 20 2018, 05:59
Цитата(ataradov @ Mar 20 2018, 03:45)
Если нужна помощь с остальным - обращайтесь, помогу.
Большое спасибо. Помощь конечно-же нужна.
Для начала немного осмыслю как все работает...
Автор: ataradov Mar 20 2018, 06:34
Какая-то из версий этой модели была на Verilog-е реализована, но исходников этого мне сейчас уже точно не найти.
Ну и не стоит эту модель воспринимать как полную истину. Это моя первая модель модема созданная по следам чтения Скляра и экспериментов.
Автор: Acvarif Mar 20 2018, 07:10
Цитата(ataradov @ Mar 20 2018, 09:34)
Какая-то из версий этой модели была на Verilog-е реализована, но исходников этого мне сейчас уже точно не найти.
Ну и не стоит эту модель воспринимать как полную истину. Это моя первая модель модема созданная по следам чтения Скляра и экспериментов.
Ну и не стоит эту модель воспринимать как полную истину. Это моя первая модель модема созданная по следам чтения Скляра и экспериментов.
Да, я понимаю.
У меня пока совсем никакого представления. Только в общих чертах. Как работает ФАПЧ мне известно. На практике приходилось это делать, но только в бинарном варианте.
Но сейчас это все реализовывается в дискретном варианте на микроконтроллерах http://electronix.ru/redirect.php?http://www.ti.com/lit/an/slaa681a/slaa681a.pdf, ПЛИС. Для начала пытаюсь хоть как-то понять это дело через Матлаб. Конечно Verilog (VHDL) код очень бы не помешал. Поскольку больше работаю с ПЛИС.
В коде появилась еще ошибка Undefined function 'diff_enc' for input arguments of type 'double'.
Похоже должна быть определена еще и функция diff_enc.
Автор: ataradov Mar 20 2018, 07:12
QUOTE (Acvarif @ Mar 20 2018, 00:10)
Похоже должна быть определена еще и функция diff_enc.
CODE
function enc=diff_enc(data, start);
if ~(exist('start'))
start = 0;
de(1)=xor(data(1), start);
for i=2:length(data), de(i)=xor(de(i-1), data(i)); end;
enc = de;
end
if ~(exist('start'))
start = 0;
de(1)=xor(data(1), start);
for i=2:length(data), de(i)=xor(de(i-1), data(i)); end;
enc = de;
end
Этот модем скорее всего перебор по сравнению с той программной реализацией на МК.
Автор: Acvarif Mar 20 2018, 07:28
Цитата(ataradov @ Mar 20 2018, 10:12)
Этот модем скорее всего перебор по сравнению с той программной реализацией на МК.
Спасибо. Дело продвинулось.
Далее Undefined function 'upsample4' for input arguments of type 'double'.
Похоже дальше есть еще неопределенные функции...
Перебор это в смысле..?
Автор: petrov Mar 20 2018, 07:30
Цитата(Acvarif @ Mar 20 2018, 10:10)
Для начала пытаюсь хоть как-то понять это дело через Матлаб. Конечно Verilog (VHDL) код очень бы не помешал. Поскольку больше работаю с ПЛИС.
Лучше сразу смотреть в строну симулинка, облегчает понимание и постепенно модель можно конретизировать до полноценного рабочего HDL кода не прибегая к внешним HDL симуляторам.
Автор: ataradov Mar 20 2018, 07:36
QUOTE (Acvarif @ Mar 20 2018, 00:28)
Похоже дальше есть еще неопределенные функции...
Вот все функции какие у меня есть.QUOTE (Acvarif @ Mar 20 2018, 00:28)
Перебор это в смысле..?
Я возможно не совсем понял статью, но они там похоже уже ожидают прямоугольный сигнал и только восстанавливают битовый клок.Этот модем еще борется с несовпадением частоты дискретизации на приемнике и передатчике (или эффектом Допплера). Плюс тут много всяких заморочек с конкретным железом. Для модема на 100 Кбод не нужна дискретизация на 60 МГц.
Автор: Acvarif Mar 20 2018, 07:51
Цитата
Лучше сразу смотреть в строну симулинка, облегчает понимание и постепенно модель можно конретизировать до полноценного рабочего HDL кода не прибегая к внешним HDL симуляторам.
Да. Согласен. Но в симуленке я сильно плаваю.
Цитата(ataradov @ Mar 20 2018, 10:36)
Вот все функции какие у меня есть.
Архив не скачивается. Выскакивает на пустую страницу.
Цитата
Я возможно не совсем понял статью, но они там похоже уже ожидают прямоугольный сигнал и только восстанавливают битовый клок.
Насколько я понял там задействован внутренний АЦП микроконтроллера. Далее цифровые библиотеки FIR и т. п.
Цитата
Этот модем еще борется с несовпадением частоты дискретизации на приемнике и передатчике (или эффектом Допплера). Плюс тут много всяких заморочек с конкретным железом. Для модема на 100 Кбод не нужна дискретизация на 60 МГц.
Доплер, это мне какраз и нужно.
Автор: ataradov Mar 20 2018, 07:55
QUOTE (Acvarif @ Mar 20 2018, 00:51)
Да. Согласен. Но в симуленке я сильно плаваю.
Может быть проще найти время и разобраться. Будет гораздо проще в итоге.QUOTE (Acvarif @ Mar 20 2018, 00:51)
Архив не скачивается. Выскакивает на пустую страницу.
Добавил снова.QUOTE (Acvarif @ Mar 20 2018, 00:51)
Насколько я понял там задействован внутренний АЦП микроконтроллера. Далее цифровые библиотеки FIR и т. п.
Одного АЦП в общем случае не достаточно нужно 2 квадратуры цифровать.Автор: Acvarif Mar 20 2018, 08:06
Цитата(ataradov @ Mar 20 2018, 10:55)
Может быть проще найти время и разобраться. Будет гораздо проще в итоге.
Согласен. Придется разбираться.
Спасибо за функции. Главный код заработал. Попытаюсь разобраться.
Цитата
Одного АЦП в общем случае не достаточно нужно 2 квадратуры цифровать.
Микроконтроллер имеет на борту многоканальный АЦП. Ну а если даже использовать один канал то при 4_х выборках на период можно использовать соседние выборки как синус и косинус.
Автор: petrov Mar 20 2018, 08:07
Цитата(Acvarif @ Mar 20 2018, 10:51)
Да. Согласен. Но в симуленке я сильно плаваю.
Получается как в анекдоте:
Мужик что-то ищет под фонарем. Тут к нему под ходит милиционер и
спрашивает: "Что вы тут делаете?" Мужик отвечает: "Ключи от квартиры
ищу". "А где потерял?". "В парке". "А зачем здесь ищешь?".
"А здесь светлее ".
Искать надо там где надо. Все эти обрывки чужих исходников - тупик. Симулинк - отличнейшее окружение, чтобы разобраться с сутью, да ещё в конце и HDL код получить можно, готовый для прошивки в FPGA.
Автор: Acvarif Mar 20 2018, 08:14
Цитата(petrov @ Mar 20 2018, 11:07)
Получается как в анекдоте:
Мужик что-то ищет под фонарем. Тут к нему под ходит милиционер и
спрашивает: "Что вы тут делаете?" Мужик отвечает: "Ключи от квартиры
ищу". "А где потерял?". "В парке". "А зачем здесь ищешь?".
"А здесь светлее ".
Искать надо там где надо. Все эти обрывки чужих исходников - тупик. Симулинк - отличнейшее окружение, чтобы разобраться с сутью, да ещё в конце и HDL код получить можно, готовый для прошивки в FPGA.
Мужик что-то ищет под фонарем. Тут к нему под ходит милиционер и
спрашивает: "Что вы тут делаете?" Мужик отвечает: "Ключи от квартиры
ищу". "А где потерял?". "В парке". "А зачем здесь ищешь?".
"А здесь светлее ".
Искать надо там где надо. Все эти обрывки чужих исходников - тупик. Симулинк - отличнейшее окружение, чтобы разобраться с сутью, да ещё в конце и HDL код получить можно, готовый для прошивки в FPGA.
Принимается.
Кстати Ваш симулинк код под названием bpsk_costas_2008_08_25.mdl надеюсь в этом поможет.
Автор: petrov Mar 20 2018, 08:29
Цитата(Acvarif @ Mar 20 2018, 11:14)
Кстати Ваш симулинк код под названием bpsk_costas_2008_08_25.mdl надеюсь в этом поможет.
Возможно, лишь как отправная точка, воспринимайте критически, ищите пробелы, недостатки и ошибки.
Автор: Acvarif Mar 22 2018, 06:26
Цитата(petrov @ Mar 20 2018, 12:29)
Возможно, лишь как отправная точка, воспринимайте критически, ищите пробелы, недостатки и ошибки.
Посмотрел модельку. Сложно для начинающего. Не понимаю, что такое скрамблер.
Google советует туториал по цифровой модуляции http://electronix.ru/redirect.php?http://www.cengage.com/resource_uploads/downloads/0495082511_315335.pdf
Только не понятно где эти все примеры находятся. В самом Матлаб, или нужно искать в сети?
Автор: ataradov Mar 22 2018, 06:32
QUOTE (Acvarif @ Mar 21 2018, 23:26)
Посмотрел модельку. Сложно для начинающего. Не понимаю, что такое скрамблер.
Нужно матлаб пока отложить и книжки по цифровой связи читать. С наскоку ничего путного не выйдет.
Обычно скрамблер убирает длинные последовательности 0 и 1 в исходном сигнале. Это нужно так как восстановление символьной частоты происходит по переходам из 0 в 1 и обратно. И если этих переходов не достаточно, то символьная синхронизация может потеряться.
Автор: Acvarif Mar 22 2018, 07:21
Цитата(ataradov @ Mar 22 2018, 10:32)
Нужно матлаб пока отложить и книжки по цифровой связи читать. С наскоку ничего путного не выйдет.
Обычно скрамблер убирает длинные последовательности 0 и 1 в исходном сигнале. Это нужно так как восстановление символьной частоты происходит по переходам из 0 в 1 и обратно. И если этих переходов не достаточно, то символьная синхронизация может потеряться.
Обычно скрамблер убирает длинные последовательности 0 и 1 в исходном сигнале. Это нужно так как восстановление символьной частоты происходит по переходам из 0 в 1 и обратно. И если этих переходов не достаточно, то символьная синхронизация может потеряться.
Спасибо. Понятно. Тогда по ходу в демодуляторе должен быть дескрамблер?
Хотя на этом думаю мне пока не стоит зацикливаться.
В простом случае принцип фазовой манипуляции ничего особенного их себя не представляет. Это все очень похоже, например, на принцип формирования кода Баркера, с чем плотно приходилось и приходится иметь дело в цифровой технике.
Что касается демодулятора с ФАПЧ (или модное название петля Костаса) это все тоже очень похоже на цифровую ФАПЧ с использованием синтезатора частоты, например, по Уолшу в обратной связи, что тоже приходилось делать, и работало оно на ура.
Думаю не стоит мне сейчас влезать в дебри цифровой связи. Это все придет по ходу.
Пока в Матлабе ищу простейший способ запустить простейший модулятор-демодулятор (кодированием или симулинк) чтобы понять как все перенести в ПЛИС. Имеется ввиду в цифровом виде (АЦП, выборки и пр.). В бинарном виде ФАПЧ ничего особенного из себя не представляет. Разве, что только пока не понимаю как подбираются промежуточная по отношению к несущей, как пересчитывать коэффициенты для управления синтезатором уолша. В свое время это все делалось на варикапах и подбиралось экспериментально, потом по мере появления новых технологий все переносилось в ПЛИС, но соблюдая все наработки.
Автор: ataradov Mar 22 2018, 07:27
QUOTE (Acvarif @ Mar 22 2018, 00:21)
Спасибо. Понятно. Тогда по ходу в демодуляторе должен быть дескрамблер?
Да, конечно.QUOTE (Acvarif @ Mar 22 2018, 00:21)
В простом случае принцип фазовой манипуляции ничего особенного их себя не представляет. Это все очень похоже, например, на принцип формирования кода Баркера, с чем плотно приходилось и приходится иметь дело в цифровой технике.
Это каша полная. Причем тут Баркер?QUOTE (Acvarif @ Mar 22 2018, 00:21)
Что касается демодулятора с ФАПЧ (или модное название петля Костаса)
Петля Костаса - это не модное название ФАПЧ, это очень специфическая схема с ФАПЧ. ФАПЧ бывает и с другими схемами, ничего общего с Костасом не имеющие.QUOTE (Acvarif @ Mar 22 2018, 00:21)
Думаю не стоит мне сейчас влезать в дебри цифровой связи. Это все придет по ходу.
Как раз таки наоборот, с пониманием цифровой связи почти автоматически придет понимание как это сделать в матлабе.Автор: Acvarif Mar 22 2018, 08:09
Цитата(ataradov @ Mar 22 2018, 11:27)
Это каша полная. Причем тут Баркер?
Баркер тоже фазовая манипуляция. Применялась для формирования сложного сигнала.
BPSK это тоже фазовая манипуляция.
Цитата
Петля Костаса - это не модное название ФАПЧ, это очень специфическая схема с ФАПЧ. ФАПЧ бывает и с другими схемами, ничего общего с Костасом не имеющие.
Не согласен. Ничего специфического там нет. В бинарном виде эта самая ФАПЧ делается с применением синтезатора частоты Уолша как гетеродина.
Используем, например, в гидроакустических лагах.
Цитата
Как раз таки наоборот, с пониманием цифровой связи почти автоматически придет понимание как это сделать в матлабе.
Согласен. Только просто чтение, по опыту, бесполезно. Пытаюсь все делать по ходу. Может и не правильно, но так уж сложилось...
Конечно хотелось-бы выяснить как все-же правильно, подобрать частоту гетеродина, например при несущей 100 кГ, если один бит при BPSK модуляции планируется передавать четырьмя периодами несущей.
Как рассчитывается шаг изменения частоты гетеродина в петле ФАПЧ.
Должны-же быть основополагающие формулы для таких расчетов...
Автор: petrov Mar 22 2018, 09:05
Цитата(Acvarif @ Mar 22 2018, 11:09)
Только просто чтение, по опыту, бесполезно. Пытаюсь все делать по ходу.
Читайте и делайте соответствующие простенькие модельки в симулинке. Про железо и предыдущий опыт сейчас надо забыть, лишнее мешает.
Цитата(Acvarif @ Mar 22 2018, 11:09)
Конечно хотелось-бы выяснить как все-же правильно, подобрать частоту гетеродина, например при несущей 100 кГ, если один бит при BPSK модуляции планируется передавать четырьмя периодами несущей.
Как рассчитывается шаг изменения частоты гетеродина в петле ФАПЧ.
Должны-же быть основополагающие формулы для таких расчетов...
Как рассчитывается шаг изменения частоты гетеродина в петле ФАПЧ.
Должны-же быть основополагающие формулы для таких расчетов...
Демодуляцию делают на нулевой частоте, частота комплексного гетеродина минус 100 кГц.
Автор: Acvarif Mar 22 2018, 11:11
Цитата(petrov @ Mar 22 2018, 13:05)
Читайте и делайте соответствующие простенькие модельки в симулинке. Про железо и предыдущий опыт сейчас надо забыть, лишнее мешает.
Пробую. Трудно дается. Неспецифическая тема.
Цитата
Демодуляцию делают на нулевой частоте, частота комплексного гетеродина минус 100 кГц.
Не совсем врубаюсь.
Комплексная это понятно.
Допустим несущая 100 кГц. Какая должна быть средняя частота управляемого синтезатора частоты DDS в петле ФАПЧ?
Тоже 100кГ?
Если взять одно квадратурное плечо-умножитель то на один вход будет подаваться модулированная несущая 100 кГц, а на другой, например sin 100 кГц. В этом случае какраз и получится нулевая чатота. Так?
Если так почему именно на нулевой, а не со сдвигом, например 100 и 102?
Мы когда-то ФАПЧ делали с небольшой разностной, например 7 и 5 кГц.
Или я чего-то не так понял?
Автор: petrov Mar 22 2018, 11:22
Цитата(Acvarif @ Mar 22 2018, 14:11)
Пробую. Трудно дается. Неспецифическая тема.
Не совсем врубаюсь.
Комплексная это понятно.
Допустим несущая 100 кГц. Какая должна быть средняя частота управляемого синтезатора частоты DDS в петле ФАПЧ?
Тоже 100кГ?
Если взять одно квадратурное плечо-умножитель то на один вход будет подаваться модулированная несущая 100 кГц, а на другой, например sin 100 кГц. В этом случае какраз и получится нулевая чатота. Так?
Если так почему именно на нулевой, а не со сдвигом, например 100 и 102?
Мы когда-то ФАПЧ делали с небольшой разностной, например 7 и 5 кГц.
Или я чего-то не так понял?
Не совсем врубаюсь.
Комплексная это понятно.
Допустим несущая 100 кГц. Какая должна быть средняя частота управляемого синтезатора частоты DDS в петле ФАПЧ?
Тоже 100кГ?
Если взять одно квадратурное плечо-умножитель то на один вход будет подаваться модулированная несущая 100 кГц, а на другой, например sin 100 кГц. В этом случае какраз и получится нулевая чатота. Так?
Если так почему именно на нулевой, а не со сдвигом, например 100 и 102?
Мы когда-то ФАПЧ делали с небольшой разностной, например 7 и 5 кГц.
Или я чего-то не так понял?
Никакие небольшие сдвиги нам не нужны. Рассматривайте перенос на нулевую, понижение частоты дискретизации отдельно от демодуляции, синхронизации, которые будут осуществляться на нулевой частоте в комплексном виде с небольшим количеством отсчётов на символ.
Автор: Acvarif Mar 22 2018, 11:46
Цитата(petrov @ Mar 22 2018, 15:22)
Никакие небольшие сдвиги нам не нужны. Рассматривайте перенос на нулевую, понижение частоты дискретизации отдельно от демодуляции, синхронизации, которые будут осуществляться на нулевой частоте в комплексном виде с небольшим количеством отсчётов на символ.
Понял. Спасибо.
Тогда если в цифровом виде (не бинарном).
Несущая 100 кгц. Частота дискретизации 400 кГц. Исходная частота NCO в петле ФАПЧ 100 кГц.
Количество периодов несущей на символ (в смысле на бит) 4.
От этого можно оттолкнуться?
Автор: petrov Mar 22 2018, 12:52
Цитата(Acvarif @ Mar 22 2018, 14:46)
Понял. Спасибо.
Тогда если в цифровом виде (не бинарном).
Несущая 100 кгц. Частота дискретизации 400 кГц. Исходная частота NCO в петле ФАПЧ 100 кГц.
Количество периодов несущей на символ (в смысле на бит) 4.
От этого можно оттолкнуться?
Тогда если в цифровом виде (не бинарном).
Несущая 100 кгц. Частота дискретизации 400 кГц. Исходная частота NCO в петле ФАПЧ 100 кГц.
Количество периодов несущей на символ (в смысле на бит) 4.
От этого можно оттолкнуться?
Нет никакой петли ФАПЧ 100 кГц, нет никакого количества периодов несущей на символ, просто переносим нужную частоту на нулевую.
Автор: Acvarif Mar 23 2018, 06:12
Цитата(petrov @ Mar 22 2018, 16:52)
Нет никакой петли ФАПЧ 100 кГц, нет никакого количества периодов несущей на символ, просто переносим нужную частоту на нулевую.
Понятно. Тогда в простейшем виде и малыми шагами.
Модулятор
|
Выход элемента Switch - BPSK сигнал
Далее, чтобы его демодулировать необходимо выделить из него комплексную огибающую. Для этого его нужно умножить на комплексную экспоненту (перенос спектра модулированного сигнала в область 0) и затем отфильтровать с помощью ФНЧ
Выглядит это так
|
Если, допустим, в модуляторе несущая равна Fn, то в демодуляторе чему должна быть равна омега0, тоже Fn?
Вот на этом моменте я почему-то застрял...
Автор: petrov Mar 23 2018, 06:49
Автор: Acvarif Mar 23 2018, 07:22
Цитата(petrov @ Mar 23 2018, 10:49)
Показан процесс демодуляции - повторное умножение на комплексную экспоненту со знаком минус.
Грубо говоря это избавление от несущей. Так?
Значит если несущая 100 кгц то в демодуляторе омега0 тоже должна быть 100 кГц. Так?
Тоесть если в квадратурном гетеродине омега0 сделать 110 кГц то спектр раздвоится на 10 и 210 кГц.
Произойдет перенос модулированного сигнала на другую частоту. А нам, грубо говоря, это н. не нужно.
Нужно всего лишь избавиться от несущей. Так?
Автор: petrov Mar 23 2018, 07:37
При умножении на комплексную экспоненту ничего не раздваивается, просто спектр сдвигается. Если мы берём реальную часть комплексного сигнала, то получаем симметричный спектр относительно нуля.
Автор: Acvarif Mar 23 2018, 08:10
Цитата(petrov @ Mar 23 2018, 11:37)
При умножении на комплексную экспоненту ничего не раздваивается, просто спектр сдвигается. Если мы берём реальную часть комплексного сигнала, то получаем симметричный спектр относительно нуля.
Да. Погорячился. Сдвигается.
Но всеравно остается вопрос почему в квадратурном демодуляторе(гетеродине) на входы умножителей обычно подается частота равная несущей?
Для того чтобы сразу от нее избавиться?
Почему не перенести спектр на другую частоту Например несущая 100 кГц Гетеродин 98 кГц. Так у нас обычно делалось. Дальше уже работа идет с 2 кГц и там уже и выделяется комплексная огибающая таким же способом перемножения комплексных експонент.
Автор: petrov Mar 23 2018, 08:17
Цитата(Acvarif @ Mar 23 2018, 11:10)
Но всеравно остается вопрос почему в квадратурном демодуляторе(гетеродине) на входы умножителей обычно подается частота равная несущей?
Вы видели созвездие BPSK, QPSK, QAM? По которым мы принимаем решения о передаваемых символах. Чтобы такую картинку получить нужно осуществить фазовую и символьную синхронизации. Если будет сдвиг по частоте, то созвездие будет вращаться с угловой частотой соответствующей сдвигу и решения принять не получится.
Автор: Acvarif Mar 23 2018, 08:51
Цитата(petrov @ Mar 23 2018, 12:17)
Вы видели созвездие BPSK, QPSK, QAM? По которым мы принимаем решения о передаваемых символах. Чтобы такую картинку получить нужно осуществить фазовую и символьную синхронизации. Если будет сдвиг по частоте, то созвездие будет вращаться с угловой частотой соответствующей сдвигу и решения принять не получится.
Спасибо. Кажется становится понятно. Нужно еще осмыслить...
Автор: Acvarif Mar 26 2018, 05:01
Для большего понимания для себя происходящего в BPSK сворганил простую модель.
CODE
%% BPSK модель
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
f= 100000;
% частота выборок
fs = f*100;
% Период выборки
Ts = 1/fs;
% Период несущей
T = 1/f;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
t = 0:Ts:n*T;
% Несущая в пакете
car = sin(2*pi*f*t);
% График данных
subplot(2, 1, 1);
stem(pData);
title('Данные');
% График несущей
subplot(2, 1, 2);
plot(car);
title('Несущая');
%% Преобразование данных в прямоугольные импульсы
tp = 0:Ts:T*M;
exdata = [];
for (i = 1:length(data))
for(j = 1:length(tp) - 1)
exdata = [exdata pData(i)];
end
end
exdata = [exdata 0];
figure;
plot(exdata, 'r-', 'LineWidth', 4);
hold on;
grid on;
plot(car, 'g-');
hold on
%% Модуляция
% Перенос exdata на несущую
mSig = exdata.*car;
plot(mSig, 'h-', 'LineWidth', 2);
hold off;
title('Прям. импульсы + несущая + промод.несущая');
%% AWGN (аддитивный белый гауссов шум) канал передачи данных
figure;
% отношение сигнал/шум в децибелах
SNR = 50;
% мощность сигнала в децибелах
SIGPOWER = 10;
% смесь мод.сигнала с шумом
rx = awgn( mSig, SNR, SIGPOWER);
plot(mSig, 'r-', 'LineWidth', 3);
hold on;
plot(rx, 'g-', 'LineWidth', 1);
grid on;
%% Демодуляция
% умножение мод.зашумленого сигн.на несущую
dem = rx.*car;
figure;
plot(dem);
grid on;
%% Декодирование
k = 1;
rcv = [];
for(i = 1:length(data))
sm = 0;
for(j = 1:length(tp) - 1)
sm = sm + dem(k);
k = k + 1;
end
if(sm > 0)
rcv = [rcv 1];
else
rcv = [rcv 0];
end
end
clc;
% Переданные данные
data
% Принятые данные
rcv
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
f= 100000;
% частота выборок
fs = f*100;
% Период выборки
Ts = 1/fs;
% Период несущей
T = 1/f;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
t = 0:Ts:n*T;
% Несущая в пакете
car = sin(2*pi*f*t);
% График данных
subplot(2, 1, 1);
stem(pData);
title('Данные');
% График несущей
subplot(2, 1, 2);
plot(car);
title('Несущая');
%% Преобразование данных в прямоугольные импульсы
tp = 0:Ts:T*M;
exdata = [];
for (i = 1:length(data))
for(j = 1:length(tp) - 1)
exdata = [exdata pData(i)];
end
end
exdata = [exdata 0];
figure;
plot(exdata, 'r-', 'LineWidth', 4);
hold on;
grid on;
plot(car, 'g-');
hold on
%% Модуляция
% Перенос exdata на несущую
mSig = exdata.*car;
plot(mSig, 'h-', 'LineWidth', 2);
hold off;
title('Прям. импульсы + несущая + промод.несущая');
%% AWGN (аддитивный белый гауссов шум) канал передачи данных
figure;
% отношение сигнал/шум в децибелах
SNR = 50;
% мощность сигнала в децибелах
SIGPOWER = 10;
% смесь мод.сигнала с шумом
rx = awgn( mSig, SNR, SIGPOWER);
plot(mSig, 'r-', 'LineWidth', 3);
hold on;
plot(rx, 'g-', 'LineWidth', 1);
grid on;
%% Демодуляция
% умножение мод.зашумленого сигн.на несущую
dem = rx.*car;
figure;
plot(dem);
grid on;
%% Декодирование
k = 1;
rcv = [];
for(i = 1:length(data))
sm = 0;
for(j = 1:length(tp) - 1)
sm = sm + dem(k);
k = k + 1;
end
if(sm > 0)
rcv = [rcv 1];
else
rcv = [rcv 0];
end
end
clc;
% Переданные данные
data
% Принятые данные
rcv
В данном случае это идеальный, полностью синхронизированный, вариант.
Я так понимаю в реальности все далеко не так. Демодулятор может иметь свою частоту дискретизации отличную от модулятора и она не синхронна с частотой дискретизации в модуляторе. Но это не самое важное. Важно что несущая подаваемая на умножитель в демодуляторе тоже не синхронна с несущей пришедшего на него сигнала, не говоря уже о доплере если источник и приемник двигаются друг относительно друга.
В этой части вопрос. Какая основная задача петли Костаса? Подстройка частоты подаваемой на умножитель, вернее на умножители демодулятора? Или задача петли Костаса это фазовая и символьная синхронизация в демодуляторе? Если это так то не врубаюсь как подстройкой частоты NCO можно добиться еще и фазовой синхронизации...
Автор: ataradov Mar 26 2018, 05:09
Костас только убирает несовпадение частот и допплер (они не различимы на приемнике). Костас ничего не делает для символьной синхронизации. Для этого в данные вставляют синхропоследовательности для начальной синхронизации и дальнейшая синхронизация осуществляется по смене фаз дынных. И чтобы этих смен было достаточно используют скремблеры.
Автор: Acvarif Mar 26 2018, 08:08
Цитата(ataradov @ Mar 26 2018, 08:09)
Костас только убирает несовпадение частот и допплер (они не различимы на приемнике). Костас ничего не делает для символьной синхронизации. Для этого в данные вставляют синхропоследовательности для начальной синхронизации и дальнейшая синхронизация осуществляется по смене фаз дынных. И чтобы этих смен было достаточно используют скремблеры.
Спасибо. Понятно. Стало более понятно значение нескольких байт преамбулы обычно используемой в пакете данных.
Попытаюсь встроить петлю Костаса в код
|
На демодулятор поступает rx (смесь модулированного сигнала с аддит. гаусовым шумом) которую далее необходимо умножать на
sin и cos поступающие с выхода VCO
Очевидно что в демодуляторе необходимо будет организовать
1. Управляемый генератор sin cos (VCO).
2. Два одинаковых ФНЧ
3. Фазовый дискриминатор
4. Еще фильтр перед управляемым генератором
Какова роль фильтра (loopfilter) после фазового дискриминатора?
Автор: petrov Mar 26 2018, 08:32
Цитата(Acvarif @ Mar 26 2018, 11:08)
Попытаюсь встроить петлю Костаса в код
В коде сложнее что-либо разглядеть, чем в симулинк модели.
Цитата(Acvarif @ Mar 26 2018, 11:08)
2. Два одинаковых ФНЧ
Согласованный фильтр лучше поставить до петли, иначе задержка фильтра будет входить в петлю ФАПЧ, с соответствующим увеличением времени настройки, для обеспечения устойчивасти.
Цитата(Acvarif @ Mar 26 2018, 11:08)
4. Еще фильтр перед управляемым генератором
Какова роль фильтра (loopfilter) после фазового дискриминатора?
Какова роль фильтра (loopfilter) после фазового дискриминатора?
Пропорционально интегрирующее звено ФАПЧ.
Ещё что-то надо делать с амплитудой сигнала, которую надо учитывать при рассчёте коэффициентов передачи ФАПЧ.
Автор: Acvarif Mar 26 2018, 10:45
Цитата(petrov @ Mar 26 2018, 11:32)
Согласованный фильтр лучше поставить до петли, иначе задержка фильтра будет входить в петлю ФАПЧ, с соответствующим увеличением времени настройки, для обеспечения устойчивасти.
До петли это как? На схеме разве не так? Фильтры стоят сразу после I Q умножителей.
Разве можно как-то иначе поставить фильтры и так чтобы не учитывалась их задерка в петле ФАПЧ?
Хотя Ваша схема ФАПЧ не похожа на привычную петлю Костаса
|
Если не сложно, поясните пожалуйста в нескольких фразах как она работает.
Цитата
Ещё что-то надо делать с амплитудой сигнала, которую надо учитывать при рассчёте коэффициентов передачи ФАПЧ.
С этим я вообще пока не знаю что делать...
В свое время когда петля ФАПЧ делалась в бинарном виде этот вопрос вообще не стоял.
Автор: petrov Mar 26 2018, 11:47
Цитата(Acvarif @ Mar 26 2018, 13:45)
До петли это как? На схеме разве не так? Фильтры стоят сразу после I Q умножителей.
Уже выше обсуждали, что комплексный сигнал у вас будет на нулевой частоте, сначала фильтр согласованный, а потом уже ФАПЧ - умножение на управляемую комплексную экспоненту, поворачиваем вращающееся созвездие в обратную сторону, чтобы оно стояло неподвижно.
Цитата(Acvarif @ Mar 26 2018, 13:45)
Хотя Ваша схема ФАПЧ не похожа на привычную петлю Костаса
Если не сложно, поясните пожалуйста в нескольких фразах как она работает.
Если не сложно, поясните пожалуйста в нескольких фразах как она работает.
Это и не петля Костаса, в оригинальной статье петля работает без символьной синхронизации, здесь же она работает после символьной синхронизации в домене одного отсчёта на символ. Согласованный фильтр здесь как раз включён в петлю.
Цитата(Acvarif @ Mar 26 2018, 13:45)
С этим я вообще пока не знаю что делать...
В свое время когда петля ФАПЧ делалась в бинарном виде этот вопрос вообще не стоял.
В свое время когда петля ФАПЧ делалась в бинарном виде этот вопрос вообще не стоял.
Моделировать, моделировать. У вас же в общем случае короткие пакеты будут приходить минимального и максимального уровня друг за другом.
Автор: Acvarif Mar 26 2018, 13:10
Цитата(petrov @ Mar 26 2018, 14:47)
Уже выше обсуждали, что комплексный сигнал у вас будет на нулевой частоте, сначала фильтр согласованный, а потом уже ФАПЧ - умножение на управляемую комплексную экспоненту, поворачиваем вращающееся созвездие в обратную сторону, чтобы оно стояло неподвижно.
Запутался.
Вот петля Костаса
|
Входной сигнал m(t) sin(ωt) Далее он сразу умножается в одном плече на sin(ωt) в другом на cos(ωt) Как результат - на выходе верхнего плеча компонента I, нижнего Q. Произошел снос на нулевую частоту. Что понимается под согласованным фильтром? У нас, например, под согласованным фильтром понимается коррелятор при работе со сложным сигналом.
Автор: petrov Mar 26 2018, 15:26
Цитата(Acvarif @ Mar 26 2018, 16:10)
Запутался.
Вот петля Костаса
Входной сигнал m(t) sin(ωt) Далее он сразу умножается в одном плече на sin(ωt) в другом на cos(ωt) Как результат - на выходе верхнего плеча компонента I, нижнего Q. Произошел снос на нулевую частоту. Что понимается под согласованным фильтром? У нас, например, под согласованным фильтром понимается коррелятор при работе со сложным сигналом.
Вот петля Костаса
Входной сигнал m(t) sin(ωt) Далее он сразу умножается в одном плече на sin(ωt) в другом на cos(ωt) Как результат - на выходе верхнего плеча компонента I, нижнего Q. Произошел снос на нулевую частоту. Что понимается под согласованным фильтром? У нас, например, под согласованным фильтром понимается коррелятор при работе со сложным сигналом.
Статья Костаса древняя, эту петлю в аналоге делали, не нужно её тупо переносить в цифру.
Разобрались, что мы можем перенести сигнал на нулевую частоту.
Далее фильтруем его согласованным фильтром, согласованным с формой видео импульсов, которые мы передаём.
Коррелятор фильтром не является, он не инвариантен ко времени прихода.
Далее петлёй Костаса на нулевой частоте разворачиваем созвездие BPSK, чтобы оно стояло.
Filter1,3 не нужны, 2 - ПИ звено. VCO - NCO, комплексная экспонента управляемая, умножается на комплексный сигнал с выхода согласованного фильтра. Детектор ошибки - произведение реальной и мнимой частей после умножителя.
Автор: Acvarif Mar 26 2018, 18:10
Спасибо что Вы мне подсказываете. Общее понимание почти сложилось. Но непонятые детали не дают покоя.
Цитата(petrov @ Mar 26 2018, 18:26)
Статья Костаса древняя, эту петлю в аналоге делали, не нужно её тупо переносить в цифру.
Да. Картинка немного не та. Вот та картинка
|
Но суть вроде не меняется.
Цитата
Разобрались, что мы можем перенести сигнал на нулевую частоту.
Железно.
Перенос на нулевую частоту это умножения входного сигнала на комплексную экспоненту.
Цитата
Далее фильтруем его согласованным фильтром, согласованным с формой видео импульсов, которые мы передаём.
В цифровом варианте это будет FIR (LPF) Это понимается под согласованным фильтром?
Цитата
Filter1,3 не нужны
Если по цифровой схеме Костаса то не нужны. Там стоят LPF.
Цитата
VCO - NCO, комплексная экспонента управляемая, умножается на комплексный сигнал с выхода согласованного фильтра.
Тут опять тупик. Получается что согласованный фильтр должен стоять до умножителей.
Может в качестве согласованного фильтра имеется ввиду это?
|
Автор: petrov Mar 26 2018, 19:15
Цитата(Acvarif @ Mar 26 2018, 21:10)
Чем отличается согласованный фильтр от ФНЧ?
Тем, что максимально возможный сигнал/шум даёт на выходе. Читайте про приподнятый косинус и корень из приподнятого косинуса.
Цитата(Acvarif @ Mar 26 2018, 21:10)
А как быть с компонентами на удвоенной несущей после переноса на нулевую?
Вроде ФНЧ их устраняет.
Вроде ФНЧ их устраняет.
Да устраняет.
Цитата(Acvarif @ Mar 26 2018, 21:10)
Тут опять тупик. Получается что
1. Согласованный фильтр должен стоять до умножителей.
2. ФНЧ не нужны.
1. Согласованный фильтр должен стоять до умножителей.
2. ФНЧ не нужны.
Да, всё уже отфильтровано, осталось только устранить вращение созвездия.
Поймите, что от переноса частоты в ноль и даунсемплинга нужно абстрагироваться, с этим вопросом можно разобраться отдельно, всё равно никакого Костаса у вас не будет в итоге, потому что ФАПЧ медленно настраивается. Далее вся обработка осуществляется на нулевой частоте, в комплексном виде, с малым количеством отсчётов на символьный интервал.
Автор: Acvarif Mar 26 2018, 19:54
Цитата(petrov @ Mar 26 2018, 22:15)
Далее вся обработка осуществляется на нулевой частоте, в комплексном виде, с малым количеством отсчётов на символьный интервал.
Да. Где-то я это понимаю. В частности про малое количество отсчетов.
Но тем не менее для моделирования должна сложиться четкая картинка компонентов и их соединения. Схема Костаса эту картинку дает.
Цитата
всё равно никакого Костаса у вас не будет в итоге, потому что ФАПЧ медленно настраивается.
Как-же тогда быть? Никакой другой схемы кроме Костаса пока не встречал...
Автор: petrov Mar 27 2018, 06:54
Цитата(Acvarif @ Mar 26 2018, 22:54)
Как-же тогда быть? Никакой другой схемы кроме Костаса пока не встречал...
Можно почитать для начала RF Architectures and Digital Signal Processing Aspects of Digital Wireless Transceivers - Nezami.
Автор: Acvarif Mar 27 2018, 07:41
Если посчитать за какое время сработает петля Костаса, например, при доплере 150 Гц (в воде это движение ~2 м/сек) (несущая 100 кГц) то получится, что для полной синхронизации потребуется время равное 150 символам (битам). Это если шаг изменения частоты NCO будет 1 Гц.
150 бит это ~18 байт преамбулы. Если длительность одного символа 4 периода несущей то это 40 мкс. Всего 150*40 = 6000 мкс. Тоесть полная подстройка произойдет за 6 мс. Это вполне нормально. Но в реальности наверняка шаг подстройки будет меньше чем 1 Гц, например 0.5 Гц. Ну тогда полная синхронизация наступит через 12 мс. Тоже вроде неплохо. Далее пойдет синхросимвол и данные.
Для случая с несущей в 100 кГц вроде неплохо.
Цитата(petrov @ Mar 27 2018, 09:54)
Можно почитать для начала RF Architectures and Digital Signal Processing Aspects of Digital Wireless Transceivers - Nezami.
Уже читаю. Глава 5.
Думаю для начала всеравно нужно построить петлю Костаса. Приступать сразу к FeedForward наверное неправильно. Ведь так?
Автор: ataradov Mar 27 2018, 07:43
VCO получает сигал ошибки в качестве контрольного, зачем перестраивать фиксированным шагом? Вам и предлагают ПИ звено поставить чтобы шаг сам нашелся какой нужно для оптимальной настройки.
Автор: Acvarif Mar 27 2018, 07:47
Цитата(ataradov @ Mar 27 2018, 10:43)
VCO получает сигал ошибки в качестве контрольного, зачем перестраивать фиксированным шагом? Вам и предлагают ПИ звено поставить чтобы шаг сам нашелся какой нужно для оптимальной настройки.
Об этом постоянно думаю. В смысле о самонастраиваемом шаге.
Что такое ПИ звено?
Автор: ataradov Mar 27 2018, 07:51
QUOTE (Acvarif @ Mar 27 2018, 00:47)
Что такое ПИ звено?
http://electronix.ru/redirect.php?https://ru.wikipedia.org/wiki/%D0%9F%D0%98%D0%94-%D1%80%D0%B5%D0%B3%D1%83%D0%BB%D1%8F%D1%82%D0%BE%D1%80
ПИ это тоже самое,только без "Д".
Автор: petrov Mar 27 2018, 08:21
Цитата(Acvarif @ Mar 27 2018, 10:41)
Уже читаю. Глава 5.
Думаю для начала всеравно нужно построить петлю Костаса.
Думаю для начала всеравно нужно построить петлю Костаса.
Это ещё месяц назад надо было сделать, в симулинке изучить все возможные воздействия и сочетания параметров, непонятно чего тянете, стройте уже. Такими темпами и 5 лет не хватит, на то, что на самом деле надо сделать для вашей задачи.
Автор: Acvarif Mar 27 2018, 12:03
Цитата(petrov @ Mar 27 2018, 11:21)
Это ещё месяц назад надо было сделать, в симулинке изучить все возможные воздействия и сочетания параметров, непонятно чего тянете, стройте уже. Такими темпами и 5 лет не хватит, на то, что на самом деле надо сделать для вашей задачи.
Ok. Попробую. В простейшем виде.
Сформировал BPSK сигнал
|
Я так понимаю сигнал необходимо оцифровать для того, чтобы потом умножать на сигнал формируемый Discrete-Time VCO.
Сделал примерно так
|
Ход верный или ерунда?
На FPGA я делал-бы именно так...
Автор: petrov Mar 27 2018, 12:19
Цитата(Acvarif @ Mar 27 2018, 15:03)
Ok. Попробую. В простейшем виде.
Сформировал BPSK сигнал
Я так понимаю сигнал необходимо оцифровать для того, чтобы потом умножать на сигнал формируемый Discrete-Time VCO.
Сделал примерно так
Ход верный или ерунда?
На FPGA я делал-бы именно так...
Сформировал BPSK сигнал
|
|
Я так понимаю сигнал необходимо оцифровать для того, чтобы потом умножать на сигнал формируемый Discrete-Time VCO.
Сделал примерно так
|
|
Ход верный или ерунда?
На FPGA я делал-бы именно так...
Не нужно моделировать сигнал на несущей, этим вы только замедляете модель, увеличивая количество отсчётов, необходимое для представления сигнала на несущей. Ни о какой оцифровке на данном этапе думать не нужно, и про FPGA тоже. Работайте с комплексными сигналами. В общем примеры я выкладывал, смотрите как там сделано.
Автор: Acvarif Mar 28 2018, 05:33
Цитата(petrov @ Mar 27 2018, 15:19)
Не нужно моделировать сигнал на несущей, этим вы только замедляете модель, увеличивая количество отсчётов, необходимое для представления сигнала на несущей. Ни о какой оцифровке на данном этапе думать не нужно, и про FPGA тоже. Работайте с комплексными сигналами. В общем примеры я выкладывал, смотрите как там сделано.
Рад-бы. Но мой уровень симулинк не позволяет сходу это сделать.
Ваш пример для меня сложен. А примеры имеющиеся в сети на мой взгляд бестолковые. Их даже в качестве рыбы нельзя использовать. Оторванные от реального железа.
Поскольку в симулинк я вообще не умею работать то он вызывает некоторое отторжение.
Петлю ФАПЧ попробую сделать в классическом матлабе.
CODE
%% BPSK модель
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
fn= 100000;
% частота выборок
fns = fn*10;
% Период выборки
Tns = 1/fns;
% Период несущей
Tn = 1/fn;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
tn = 0:Tns:n*Tn;
% Несущая в пакете
car = sin(2*pi*fn*tn);
%% Преобразование данных в прямоугольные импульсы
tpn = 0:Tns:Tn*M;
exdata = [];
for (i = 1:length(data))
for(j = 1:length(tpn) - 1)
exdata = [exdata pData(i)];
end
end
exdata = [exdata 0];
%% Модуляция
% Перенос exdata на несущую
mSig = exdata.*car;
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
fn= 100000;
% частота выборок
fns = fn*10;
% Период выборки
Tns = 1/fns;
% Период несущей
Tn = 1/fn;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
tn = 0:Tns:n*Tn;
% Несущая в пакете
car = sin(2*pi*fn*tn);
%% Преобразование данных в прямоугольные импульсы
tpn = 0:Tns:Tn*M;
exdata = [];
for (i = 1:length(data))
for(j = 1:length(tpn) - 1)
exdata = [exdata pData(i)];
end
end
exdata = [exdata 0];
%% Модуляция
% Перенос exdata на несущую
mSig = exdata.*car;
Теория говорит, что mSig по сути комплексный сигнал. Но в демодуляторе нужно использовать только реальную его часть, которую затем в петле ФАПЧ нужно умножать на комплексную экспоненту.
Какой выглядит в коде реальная часть mSig?
Автор: ataradov Mar 28 2018, 06:47
Нужно взять исходный бинарный сигнал, умножить его на синус и косинус частоты Допплера без всяких несущих. Все, получили I и Q. Это сигнал идущий на схему вращающую созвездие. Костас или не Костас - это как захочется.
Чтобы визуально увидеть созвездие нужно построить график I(Q). Для нулевого Допплера и в отсутствии шума и прочих искажений это будут 2 точки. При добавлении Допплера это будет окружность. Задача ФАПЧ вернуть эти точки на исходное место и следить, чтобы они там и остались.
Автор: Acvarif Mar 28 2018, 07:09
Цитата(ataradov @ Mar 28 2018, 09:47)
Нужно взять исходный бинарный сигнал, умножить его на синус и косинус частоты Допплера без всяких несущих. Все, получили I и Q. Это сигнал идущий на схему вращающую созвездие. Костас или не Костас - это как захочется.
Чтобы визуально увидеть созвездие нужно построить график I(Q). Для нулевого Допплера и в отсутствии шума и прочих искажений это будут 2 точки. При добавлении Допплера это будет окружность. Задача ФАПЧ вернуть эти точки на исходное место и следить, чтобы они там и остались.
Чтобы визуально увидеть созвездие нужно построить график I(Q). Для нулевого Допплера и в отсутствии шума и прочих искажений это будут 2 точки. При добавлении Допплера это будет окружность. Задача ФАПЧ вернуть эти точки на исходное место и следить, чтобы они там и остались.
Понял. Спасибо.
Автор: Acvarif Mar 30 2018, 05:13
Не получается.
Пытаюсь сформировать синфазную и квадратурную составляющую модулирующего сигнала.
CODE
%% BPSK модель
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
fn = 100000;
% Доплер
fd = 100;
% частота выборок
fns = fn*8;
% Период выборки
Tns = 1/fns;
% Период несущей
Tn = 1/fn;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
tn = 0:Tns:n*Tn;
% Несущая в пакете
car_sin = sin(2*pi*fn*tn);
car_cos = cos(2*pi*fn*tn);
%% Преобразование данных в прямоугольные импульсы
tpn_i = 0:Tns:Tn*M;
exdata_i = [];
for (i = 1:length(data))
for(j = 1:length(tpn_i) - 1)
exdata_i = [exdata_i pData(i)];
end
end
exdata_i = [exdata_i 0];
tpn_q = 0:Tns:Tn*M;
exdata_q = [];
for (i = 1:length(data))
for(j = 1:length(tpn_q) - 1)
% exdata_q = [exdata_q pData(i)];
if (length(tpn_q) == 17)
exdata_q = [exdata_q 1];
else
exdata_q = [exdata_q -1];
end
end
end
exdata_q = [exdata_q 0];
%% Модуляция
% Перенос exdata на несущую
mSig_i = exdata_i.*car_cos;
mSig_q = exdata_q.*car_sin;
%%
figure;
plot(exdata_i, 'r-', 'LineWidth', 4);
hold on;
grid on;
plot(exdata_q, 'b-', 'LineWidth', 4);
hold on;
plot(car_sin, 'g-');
hold on
hold off;
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
fn = 100000;
% Доплер
fd = 100;
% частота выборок
fns = fn*8;
% Период выборки
Tns = 1/fns;
% Период несущей
Tn = 1/fn;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
tn = 0:Tns:n*Tn;
% Несущая в пакете
car_sin = sin(2*pi*fn*tn);
car_cos = cos(2*pi*fn*tn);
%% Преобразование данных в прямоугольные импульсы
tpn_i = 0:Tns:Tn*M;
exdata_i = [];
for (i = 1:length(data))
for(j = 1:length(tpn_i) - 1)
exdata_i = [exdata_i pData(i)];
end
end
exdata_i = [exdata_i 0];
tpn_q = 0:Tns:Tn*M;
exdata_q = [];
for (i = 1:length(data))
for(j = 1:length(tpn_q) - 1)
% exdata_q = [exdata_q pData(i)];
if (length(tpn_q) == 17)
exdata_q = [exdata_q 1];
else
exdata_q = [exdata_q -1];
end
end
end
exdata_q = [exdata_q 0];
%% Модуляция
% Перенос exdata на несущую
mSig_i = exdata_i.*car_cos;
mSig_q = exdata_q.*car_sin;
%%
figure;
plot(exdata_i, 'r-', 'LineWidth', 4);
hold on;
grid on;
plot(exdata_q, 'b-', 'LineWidth', 4);
hold on;
plot(car_sin, 'g-');
hold on
hold off;
Получается ерунда.
У Вас похоже это работает
CODE
clear all;
Br = 117.1875*10^3;
Fd = Br*8;
Fif = 20*10^6;
Fdif = 60*10^6;
Fdopl = 100;
%% === TRANSMITTER ===
M = pngen(5);
%M = 0.5;
D = rand(1, 200*8) > 0.5;
%D = [1 0 1 0];
CRC = zeros(1, 16);
for i=1:length(D),
CRC(16) = xor(CRC(16), D(i));
CRC = [CRC(2:16) CRC(1)];
end;
FILL = ones(1, 200)*0.5;
DATA = [dec2bin(length(D)/8, 11)=='1' D CRC];
packet = [FILL M M diff_enc([0 DATA]) FILL]*2-1;
% redescr
Ndkr = Fd/Br;
pkt = packet(ceil(0.01:1/Ndkr:length(packet)));
hch = firls(31, [0 .24 .3 .5]*2, [1 1 0 0]);
pktf = filter(hch, [1], pkt);
% upsample
t1 = upsample4(pktf, 31);
t2 = upsample4(t1, 31);
pktfu = upsample4(t2, 31)';
% -> IF
l = length(pktfu);
s = sin(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
c = cos(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
pkts = pktfu.*s(1:l);
pktc = pktfu.*c(1:l);
%% === TRANSMIT LINE ===
Ir = awgn(awgn(pkts, 0.05), 0.01);
Qr = awgn(awgn(pktc, 0.05), 0.01);
%% === GRAPHICS ===
figure;
plot(pkts, 'g-');
grid on;
hold on
plot(pktc, 'r-');
hold on
figure;
plot(Ir, 'g-');
grid on;
hold on
plot(Qr, 'r-');
hold on
plot(pktfu, 'r-');
hold on
figure;
plot(pkt, 'g-');
grid on;
hold on
clc;
% Переданные данные
%pkt
Br = 117.1875*10^3;
Fd = Br*8;
Fif = 20*10^6;
Fdif = 60*10^6;
Fdopl = 100;
%% === TRANSMITTER ===
M = pngen(5);
%M = 0.5;
D = rand(1, 200*8) > 0.5;
%D = [1 0 1 0];
CRC = zeros(1, 16);
for i=1:length(D),
CRC(16) = xor(CRC(16), D(i));
CRC = [CRC(2:16) CRC(1)];
end;
FILL = ones(1, 200)*0.5;
DATA = [dec2bin(length(D)/8, 11)=='1' D CRC];
packet = [FILL M M diff_enc([0 DATA]) FILL]*2-1;
% redescr
Ndkr = Fd/Br;
pkt = packet(ceil(0.01:1/Ndkr:length(packet)));
hch = firls(31, [0 .24 .3 .5]*2, [1 1 0 0]);
pktf = filter(hch, [1], pkt);
% upsample
t1 = upsample4(pktf, 31);
t2 = upsample4(t1, 31);
pktfu = upsample4(t2, 31)';
% -> IF
l = length(pktfu);
s = sin(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
c = cos(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
pkts = pktfu.*s(1:l);
pktc = pktfu.*c(1:l);
%% === TRANSMIT LINE ===
Ir = awgn(awgn(pkts, 0.05), 0.01);
Qr = awgn(awgn(pktc, 0.05), 0.01);
%% === GRAPHICS ===
figure;
plot(pkts, 'g-');
grid on;
hold on
plot(pktc, 'r-');
hold on
figure;
plot(Ir, 'g-');
grid on;
hold on
plot(Qr, 'r-');
hold on
plot(pktfu, 'r-');
hold on
figure;
plot(pkt, 'g-');
grid on;
hold on
clc;
% Переданные данные
%pkt
Но там слишком большое количество данных. Сложно рассмотреть подробности.
Автор: ataradov Mar 30 2018, 05:28
Для отладки синхронизации не нужна несущая, ее можно добавить в самом конце.
Вот как будут выглядеть I и Q компоненты идеального сигнала после всех переносов частоты.
CODE
clear all;
Br = 400e3;
Fd = Br*8;
Fdopl = 100;
D = (rand(1, 10000) > 0.5) * 2 - 1;
% redescr
Ndkr = Fd/Br;
pkt = D(ceil(0.01:1/Ndkr:length(D)));
l = length(pkt);
s = sin(2*pi*Fdopl*(0:1/Fd:l/Fd));
c = cos(2*pi*Fdopl*(0:1/Fd:l/Fd));
I = s(1:l) .* pkt;
Q = c(1:l) .* pkt;
figure; hold on; plot(I); plot(Q, 'r');
figure; hold on; plot(I, Q, '*');
Br = 400e3;
Fd = Br*8;
Fdopl = 100;
D = (rand(1, 10000) > 0.5) * 2 - 1;
% redescr
Ndkr = Fd/Br;
pkt = D(ceil(0.01:1/Ndkr:length(D)));
l = length(pkt);
s = sin(2*pi*Fdopl*(0:1/Fd:l/Fd));
c = cos(2*pi*Fdopl*(0:1/Fd:l/Fd));
I = s(1:l) .* pkt;
Q = c(1:l) .* pkt;
figure; hold on; plot(I); plot(Q, 'r');
figure; hold on; plot(I, Q, '*');
Начинайте с этого кода и добавляйте схемы, которые восстановят этот сигнал.
Если запустить с Fdopl = 0, то на второй картинке будет 2 звездочки. Это идеальный конечный результат. Если запустить Fdopl = 100, то будет окружность, образованная вращением этих звездочек.
Задача ФАПЧ вернуть звездочки на место и поддерживать их там.
Как только заработает с идеальным сигналом, то нужно будет обрезать ему полосу и добавить шум. И отлаживать это. И только после этого можно смотреть на несущие, если сильно хочется.
Автор: Acvarif Mar 30 2018, 06:51
Цитата(ataradov @ Mar 30 2018, 08:28)
Для отладки синхронизации не нужна несущая, ее можно добавить в самом конце.
Вот как будут выглядеть I и Q компоненты идеального сигнала после всех переносов частоты.
Вот как будут выглядеть I и Q компоненты идеального сигнала после всех переносов частоты.
Спасибо. Попробую разобраться.
Автор: Acvarif Apr 12 2018, 10:44
Сворганил BPSK приемо-передатчик с петлей Костаса на базе VCO.
Основной код петли заимствовал на просторах сети.
CODE
%% Петля Костаса на базе VCO
%% Модуляция
clear;
% Битовая картика профиля pacman
image = [ 1, 1, 1, 0, 0, 1, 1, 1;
1, 1, 0, 0, 0, 0, 1, 1;
1, 0, 0, 0, 1, 0, 0, 1;
0, 0, 0, 0, 0, 0, 1, 1;
0, 0, 0, 0, 0, 1, 1, 1;
1, 0, 0, 0, 0, 0, 0, 1;
1, 1, 0, 0, 0, 0, 1, 1;
1, 1, 1, 0, 0, 1, 1, 1 ];
% Начальный символ приема
sync_symbol = [ 1 ];
% Битовый поток с начальным символом
bit_stream = [ sync_symbol ];
% Генерация битового потока из матрицы профиля pacman
for j = 1:8
bit_stream = [bit_stream image(j,
];end
%bit_stream = (rand(1, 6) > 0.5);
%bit_stream = [1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0];
bit_stream = [1 1 0 1 0];
N = length(bit_stream);
fprintf('bit_stream %d\n', bit_stream);
% Частота семплирования 288 кГц
fs = 288000;
% Несущая частота относительно частоты выборок 36 кГц
f = 0.125;
% Частота доплера 300 Гц
fdop = f/1200;
% Символьная частота 100 Гц
fsim = 100;
% Количество выборок на период несущей 8
ns = 1/f;
% Количество выборок необх. для одного символа (360*8 периодов нес. f)
% Символьная скорость 100 Гц
tsymbol = fs*f*ns/fsim;
fprintf('tsymbol %d\n', tsymbol);
% Счетчик выборок в 0
t = 0;
% Количество символов переданных с нулевой фазой
nDummySymbols = 4;
% Общее количество выборок
%ttotal = N*tsymbol + nDummySymbols * tsymbol + tsymbol*3;
ttotal = N*tsymbol + nDummySymbols * tsymbol;
fprintf('ttotal %d\n', ttotal);
% Создать пустой сигнал
BPSK_signal = zeros(ttotal,1);
% Символы с нулевой фазой
for i = 1:nDummySymbols
for j = 1:tsymbol
BPSK_signal(t+1) = cos( 2*pi*(f+fdop)*t );
BPSK_signal_nomod(t+1) = cos(2*pi*(f+fdop)*t);
t = t + 1;
end
end
fprintf('t %d\n', t);
% BPSK сигнал
for i = 1:length(bit_stream)
phase = pi * bit_stream(i);
for j = 1:tsymbol
BPSK_signal(t+1) = cos(2*pi*(f+fdop)*t + phase);
BPSK_signal_mod(t+1) = cos(2*pi*(f+fdop)*t + phase);
BPSK_signal_nomod(t+1) = 0;
t = t + 1;
end
end
fprintf('length BPSK_signal %d\n', length(BPSK_signal));
fprintf('t %d\n', t);
%% AWGN (аддитивный белый гауссов шум) канал передачи данных
% отношение сигнал/шум в децибелах
SNR = 20;
% мощность сигнала в децибелах
SIGPOWER = 10;
% смесь мод.сигнала с шумом
BPSK_signal = awgn( BPSK_signal, SNR, SIGPOWER);
%% Демодуляция
% полученная картинка в 0
rec_image = zeros(8,8);
% Координаты картинки
x = 1;
y = 1;
% начальный шаг в 0
t = 0;
% Все матрицы в 0
carrier_inph = zeros(ttotal,1);
mixer_inph = zeros(ttotal,1);
demod_thr_inph = zeros(ttotal,1);
% Фильтр демодулятора
b = fir1(100, 0.1);
% Память состояния демодулятора
stateLowpassInph = zeros(numel(
-1,1);% переменные для PLL осциллятора
% cosine output
c = 1;
% delayed cosine by one timestep
c_delay = 0;
% sine output
s = 0;
% delayed sine by one timestep
s_delay = 0;
% Данные из VCO для целей отладки
sine = zeros(ttotal,1);
cosine =zeros(ttotal,1);
vco = zeros(ttotal,1);
% PLL loop фильтр
bPLL = fir1(40, 0.1);
% and the state memory
zfPLLa = zeros(numel(bPLL)-1,1);
zfPLLc = zeros(numel(bPLL)-1,1);
zfPLLs = zeros(numel(bPLL)-1,1);
% Обратный отсчет. При нулевом значении отбирается символ.
% Начальный обратный отсчет учитывает задержку
% fir lowpass и пропускает начальный символ
sampleMoment = tsymbol/2+tsymbol;
% Это время необходимое PLL для стабилизации на несущей
nPLLsettle = (nDummySymbols/2) * tsymbol;
% Игнорировать данные, пока не получен символ "1"
ignoreData = 1;
% VCO шаг
VCOgain = 0.005;
% VCO, которое управляет частотой. при v = 0 он точно равен f.
v = 0;
% Основной цикл
for step = 1:ttotal
% Получить значение BPSK_signal
sample = BPSK_signal(step);
% PLL - "voltage" управляемый осциллятор с +/- v входом
% w0 частота VCO
w0 = 2*pi*(f+v*VCOgain);
% Cохранить предыдущее состояние
c_delay = c;
s_delay = s;
% Подсчет нового состояния
c = c_delay * cos(w0) - s_delay * sin(w0);
s = s_delay * cos(w0) + c_delay * sin(w0);
% Cохраним все в удобных векторах для графики
sine(step) = s;
cosine(step) = c;
vco(step) = v;
% конец VCO
% Сохраняем sine/cosine. Обратите внимание, что инфазная волна
% теперь является синусоидальной волной, потому что PLL с мультипликацией
% в качестве фазового детектора имеет сдвиг фазы на 90 градусов между
% входом и выходом
carrier_inph(step) = c;
% Это решает, находимся ли мы в режиме синхронизации PLL или режима замораживания
% Примечание: это должно быть заменено блокировкой PLL в реальном приеме,
% потому что мы не знаем, когда передача фактически
% начинается. Здесь мы просто предполагаем, что она начинается с самого начала
% PLL in action: изменение частоты через v
% Фазовый детектор
pc = sample*c;
ps = sample*-s;
% filtering out the 2f term
% Фильтрация pc
[vc,zfPLLc] = filter(bPLL,1,pc,zfPLLc);
% filtering out the 2f term
% Фильтрация ps
[vs,zfPLLs] = filter(bPLL,1,ps,zfPLLs);
[v,zfPLLa] = filter(bPLL,1,vc*vs,zfPLLa);
% we sneak here a -1 in because this is ambigous
% шаг на -1
mixer_inph(step) = -sample*carrier_inph(step);
[lp_demod_inph(step),stateLowpassInph] = filter(b,1,mixer_inph(step),stateLowpassInph);
demod_thr_inph(step) = lp_demod_inph(step);
if (ignoreData)
% Ждем стартового символа 1
if ( lp_demod_inph(step) > 0.25 )
% got it!
% Стартовый символ прошел
ignoreData = 0;
end
else
% Ожидание середины символа для выборки демодулятора
sampleMoment = sampleMoment - 1;
if (sampleMoment == 0)
sampleMoment = tsymbol;
if ( demod_thr_inph(step) > 0.2 )
lsb = 1;
else
lsb = 0;
end
fprintf('lsb %d\n', lsb);
% debug samples: add small spikes
% отладка: добавка небольших пиков
lp_demod_inph(step) = lp_demod_inph(step) + 0.1;
% Сохранить результат в картинку
if (((x<9) && (y<9)) && (ignoreData == 0))
rec_image(y,x) = lsb;
x = x + 1;
if (x > 8)
x = 1;
y = y + 1;
end % if (x > 8)
end % x>9 && y>9
end % sampleMoment == 0
end % ignore data
% Следующий шаг в выборке
t = t + 1;
end
%% Графика
% Немодулированный сигнал
figure;
plot(BPSK_signal_nomod, 'b-', 'LineWidth', 2);
hold on;
grid on;
title('Немодулированный сигнал');
% Модулированный сигнал
figure;
plot(BPSK_signal_mod, 'b-', 'LineWidth', 2);
hold on;
grid on;
title('Модулированный сигнал');
% Полный BPSK_signal
figure;
plot(BPSK_signal, 'b-', 'LineWidth', 2);
hold on;
grid on;
title('BPSK сигнал');
% Демодулированный сигнал
figure;
subplot(2,1,1);
plot(lp_demod_inph,'LineWidth',2);
xlabel('Samples');
ylabel('Amplitude');
title('demod / lowpass filtered');
grid on;
% Картинка
subplot(2,1,2);
imshow(rec_image);
% Напряжение VCO
figure;
subplot(1,1,1);
plot(vco);
Работает при соотнощении сигнал-шум 20 на 10 децибел, доплере 30 Гц при несущей в 36 кГц и символьной скорости 100 Гц
Напрягает то, что при соотношении сигнал шум 1 приемник перестает работать. Все разваливается, даже при нулевом доплере.
Все так и должно быть или нет?
Получается что для нормальной работы BPSK приемника необходим хороший приемный усилитель с узкой полосой (определяется удвоенной символьной скоростью) и коэффициентом усиления достаточным, чтобы при минимально возможном сигнале на его входе он мог обеспечить необходимый уровень сигнала на входе BPSK приемника.
Но с другой стороны при максимально возможном уровне сигнала на входе приемного усилителя он не должен входить в ограничение выхода.
Не представляю как решается такая задачка.
Автор: Acvarif Apr 16 2018, 07:54
Ниже код BPSK модулятор-демодулятор с петлей Костаса на базе VCO для акустического модема.
Просьба покритиковать. На базе данного кода предполагается писать код для ПЛИС либо микроконтроллера, только с использованием NCO.
Характеристики:
- Символьная скорость 100 бит в сек
- Несущая частота 36 000 Гц
- 8 выборок на период несущей
- Устойчиво работает при соотношении сигнал-шум 1
- Максимальный Доплер 60 Гц (в воде ~ 2.5 м/сек)
Есть сомнения в необходимости входного полосового фильтра поскольку наверняка некоторую селективность будет иметь аналоговый приемный тракт.
CODE
%% Петля Костаса на базе VCO
%% Модуляция
clear;
% Битовая картика профиля pacman
% image = [ 1, 1, 1, 0, 0, 1, 1, 1;
% 1, 1, 0, 0, 0, 0, 1, 1;
% 1, 0, 0, 0, 1, 0, 0, 1;
% 0, 0, 0, 0, 0, 0, 1, 1;
% 0, 0, 0, 0, 0, 1, 1, 1;
% 1, 0, 0, 0, 0, 0, 0, 1;
% 1, 1, 0, 0, 0, 0, 1, 1;
% 1, 1, 1, 0, 0, 1, 1, 1 ];
image = [ 1, 0, 1, 0, 1, 0, 1, 0;
0, 1, 0, 1, 0, 1, 0, 1;
1, 0, 1, 0, 1, 0, 1, 0;
0, 1, 0, 1, 0, 1, 0, 1;
1, 1, 1, 1, 1, 1, 1, 1;
0, 0, 0, 0, 0, 0, 0, 0;
1, 1, 1, 1, 1, 1, 1, 1;
0, 0, 0, 0, 0, 0, 0, 0 ];
% Начальный символ приема
sync_symbol = [ 1 ];
% Битовый поток с начальным символом
bit_stream = [ sync_symbol ];
% Генерация битового потока из матрицы профиля pacman
for j = 1:8
bit_stream = [bit_stream image(j,];
end
%bit_stream = (rand(1, 20) > 0.5);
%bit_stream = [1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0];
%bit_stream = [1 1 0 1 0 1 0 1 0];
N = length(bit_stream);
fprintf('bit_stream %d\n', bit_stream);
% Частота семплирования 288 кГц (8 выборок на период)
fs = 288000;
% Несущая частота относительно частоты выборок 36 кГц
f = 0.125;
% Частота доплера 60 Гц (~2.5 м.сек)
fdop = f/600;
% % Частота семплирования 144 кГц (4 выборки на период)
% fs = 144000;
% % Несущая частота относительно частоты выборок 36 кГц
% f = 0.25;
% % Частота доплера 60 Гц (~2.5 м/сек)
% fdop = -f/600;
% Частота среза НЧ фильтра в демодуляторе
fc = 100;
% Частота Найквиста
fn = fs/2;
% Порядок НЧ фильтра демодулятора
lowpfiltord = 50;
% Порядок loop фильтра демодулятора
loopfiltord = 50;
% Амплитуда
A = 1.0;
% Символьная скорость 100 Гц
fsim = 100;
% Количество выборок на период несущей
ns = 1/f;
% Количество выборок необх. для одного символа (360*8 периодов нес. f)
% при символьной скорости fsim
tsymbol = fs*f*ns/fsim;
fprintf('tsymbol %d\n', tsymbol);
% Счетчик выборок в 0
t = 0;
% Количество символов переданных с нулевой фазой (пилот)
nDummySymbols = 4;
% Общее количество выборок
%ttotal = N*tsymbol + nDummySymbols * tsymbol + tsymbol*3;
ttotal = N*tsymbol + nDummySymbols * tsymbol;
fprintf('ttotal %d\n', ttotal);
% Создать пустой BPSK сигнал
BPSK_signal_t = zeros(ttotal,1);
% Пилот сигнал
for i = 1:nDummySymbols
for j = 1:tsymbol
BPSK_signal_t(t+1) = A*sin( 2*pi*(f+fdop)*t );
BPSK_signal_nomod(t+1) = A*sin(2*pi*(f+fdop)*t);
t = t + 1;
end
end
fprintf('t %d\n', t);
% BPSK сигнал
for i = 1:length(bit_stream)
phase = pi * bit_stream(i);
for j = 1:tsymbol
BPSK_signal_t(t+1) = A*sin(2*pi*(f+fdop)*t + phase);
BPSK_signal_mod(t+1) = A*sin(2*pi*(f+fdop)*t + phase);
BPSK_signal_nomod(t+1) = 0;
t = t + 1;
end
end
fprintf('length BPSK_signal %d\n', length(BPSK_signal_t));
fprintf('t %d\n', t);
%% AWGN (аддитивный белый гауссов шум) + полосовая фильтрация
% отношение сигнал/шум в децибелах
SNR = 9;
% мощность сигнала в децибелах
SIGPOWER = 12;
% смесь мод.сигнала с шумом
BPSK_signal_awgn = awgn( BPSK_signal_t, SNR, SIGPOWER );
% Порядок полосового фильтра
bandpfiltord = 50;
% нижняя частота полосового фильтра 35КГц
W1=35e3/(0.5*fs);
% верхняя частота полосового фильтра 37КГц
W2=37e3/(0.5*fs);
% расчет окна Ханна
k=hann(bandpfiltord+1);
% расчет коэффициентов входного полосового фильтра
bp=fir1(bandpfiltord, [W1 W2], k);
% Память состояния фильтра
stateBandpassBpsk = zeros(numel(bp)-1,1);
% Полосовая фильтрация BPSK сигнала + AWGN шум
[BPSK_signal_awgn_filtr,stateBandpassBpsk] = filter(bp,1,BPSK_signal_awgn,stateBandpassBpsk);
% BPSK_signal_awgn_filtr = filter(bp, 1, BPSK_signal_awgn);
% BPSK_signal
%BPSK_signal = BPSK_signal_awgn;
BPSK_signal = BPSK_signal_awgn_filtr;
%% Демодуляция
% полученная картинка в 0
rec_image = zeros(8,8);
% Координаты картинки
x = 1;
y = 1;
% начальный шаг в 0
t = 0;
% Все матрицы в 0
carrier_inph = zeros(ttotal,1);
mixer_inph = zeros(ttotal,1);
demod_thr_inph = zeros(ttotal,1);
%расчет окна Кайзера
kb=kaiser(lowpfiltord+1);
% Коэффициенты НЧ фильтра демодулятора порядок lowpfiltord
b = fir1(lowpfiltord, fc/fn, kb);
% Память состояния для НЧ фильтра демодул.
stateLowpassInph = zeros(numel(-1,1);
% переменные для PLL осциллятора
% cosine output
c = 1;
% delayed cosine by one timestep
c_delay = 0;
% sine output
s = 0;
% delayed sine by one timestep
s_delay = 0;
% Данные из VCO для целей отладки
sine = zeros(ttotal,1);
cosine =zeros(ttotal,1);
vco = zeros(ttotal,1);
%расчет окна Ханна
kp=hann(loopfiltord+1);
% Коэффициенты PLL loop фильтра
bPLL = fir1(loopfiltord, fc/fn, kp);
% Обнулить память
zfPLLa = zeros(numel(bPLL)-1,1);
zfPLLc = zeros(numel(bPLL)-1,1);
zfPLLs = zeros(numel(bPLL)-1,1);
% Обратный отсчет. При нулевом значении отбирается символ.
% Начальный обратный отсчет учитывает задержку
% fir lowpass и пропускает начальный символ
sampleMoment = tsymbol/2+tsymbol;
% Это время необходимое PLL для стабилизации на несущей
nPLLsettle = (nDummySymbols/2) * tsymbol;
% Игнорировать данные, пока не получен символ "1" (синхробит)
ignoreData = 1;
% VCO шаг
VCOgain = 0.005;
% VCO, которое управляет частотой. при v = 0 он точно равен f.
v = 0;
% Основной цикл
for step = 1:ttotal
% Полосовая фильтрация BPSK сигнала + AWGN шум
% [BPSK_signal(step),stateBandpassBpsk] = filter(bp,1,BPSK_signal_awgn(step),stateBandpassBpsk);
% Получить значение BPSK_signal в sample
sample = BPSK_signal(step);
% PLL - управляемый осциллятор с +/- v входом
% w0 частота VCO
w0 = 2*pi*(f+v*VCOgain);
% Cохранить предыдущее состояние
c_delay = c;
s_delay = s;
% Подсчет нового состояния
c = c_delay * cos(w0) - s_delay * sin(w0);
s = s_delay * cos(w0) + c_delay * sin(w0);
% Cохраним все в удобных векторах для графики
sine(step) = s;
cosine(step) = c;
vco(step) = v;
% конец VCO
% Сохраняем sine/cosine. Обратите внимание, что инфазная волна
% теперь является синусоидальной волной, потому что PLL с мультипликацией
% в качестве фазового детектора имеет сдвиг фазы на 90 градусов между
% входом и выходом
carrier_inph(step) = c;
% Это решает, находимся ли мы в режиме синхронизации PLL или режиме замораживания
% Примечание: это должно быть заменено блокировкой PLL в реальном приеме,
% потому что мы не знаем, когда передача фактически
% начинается. Здесь мы просто предполагаем, что она начинается с самого начала
% PLL изменение частоты через v
% Фазовый детектор
pc = sample*c;
ps = sample*-s;
% Фильтрация pc
[vc,zfPLLc] = filter(bPLL,1,pc,zfPLLc);
% Фильтрация ps
[vs,zfPLLs] = filter(bPLL,1,ps,zfPLLs);
[v,zfPLLa] = filter(bPLL,1,vc*vs,zfPLLa);
% шаг на -1
mixer_inph(step) = sample*carrier_inph(step);
[lp_demod_inph(step),stateLowpassInph] = filter(b,1,mixer_inph(step),stateLowpassInph);
demod_thr_inph(step) = lp_demod_inph(step);
if (ignoreData)
% Ждем стартового символа 1
if ( lp_demod_inph(step) > 0.25 )
% Стартовый символ прошел
ignoreData = 0;
end
else
% Ожидание середины символа для выборки демодулятора
sampleMoment = sampleMoment - 1;
if (sampleMoment == 0)
sampleMoment = tsymbol;
if ( demod_thr_inph(step) > 0.1 )
lsb = 1;
else
lsb = 0;
end
fprintf('lsb %d\n', lsb);
% отладка: добавка небольших пиков
lp_demod_inph(step) = lp_demod_inph(step) + 0.1;
% Сохранить результат в картинку 8:8
if (((x<9) && (y<9)) && (ignoreData == 0))
rec_image(y,x) = lsb;
x = x + 1;
if (x > 8)
x = 1;
y = y + 1;
end % if (x > 8)
end % x>9 && y>9
end % sampleMoment == 0
end % ignoreData
% Следующий шаг в выборке
t = t + 1;
end
%% Графика
% % Немодулированный сигнал
% figure;
% subplot(2,1,1);
% plot(BPSK_signal_nomod, 'b-', 'LineWidth', 2);
% hold on;
% grid on;
% title('Немодулированный сигнал (пилот)');
%
% % Модулированный сигнал
% %figure;
% subplot(2,1,2);
% plot(BPSK_signal_mod, 'r-', 'LineWidth', 2);
% hold on;
% grid on;
% title('Модулированный сигнал (данные)');
% Полный BPSK_signal
figure;
subplot(2,1,1);
plot(BPSK_signal_awgn, 'b-', 'LineWidth', 2);
xlabel('Выборки');
hold on;
grid on;
title('Полный BPSK сигнал + шум AWGN');
subplot(2,1,2);
plot(BPSK_signal, 'r-', 'LineWidth', 2);
xlabel('Выборки');
hold on;
grid on;
title('Полный BPSK сигнал + шум AWGN + Фильтр');
% Демодулированный сигнал
figure;
subplot(2,1,1);
plot(lp_demod_inph,'LineWidth',2);
xlabel('Выборки');
ylabel('Амплитуда');
title('Демодулированный сигнал / НЧ фильтрация');
grid on;
% Картинка
subplot(2,1,2);
imshow(rec_image);
title('Битовая картинка 8:8 квадратиков 1-белый, 0-черный кв. ' );
% Напряжение VCO
% figure
% plot(vco);
% title('Напряжение VCO');
%
% % Фильтр демодулятора
% figure
% % АЧХ НЧ фильтра
% [h,f]=freqz(b,1);
% % Параметры являются соответственно частотой и амплитудой
% plot(f*fs/(2*pi),20*log10(abs(h)))
% xlabel('частота/Hz');ylabel('амплитуда/dB');title('АЧХ НЧ фильтра демодулятора');
%
% % Фильтр полосовой на 36 кГц
% figure
% [h,fp]=freqz(bp,1);%amplitude-frequency characteristic graph
% plot(fp*fs/(2*pi),20*log10(abs(h)))%parameters are respectively frequency and amplitude
% xlabel('частота/Hz');ylabel('амплитуда/dB');title('АЧХ полосового фильтра');
%% Модуляция
clear;
% Битовая картика профиля pacman
% image = [ 1, 1, 1, 0, 0, 1, 1, 1;
% 1, 1, 0, 0, 0, 0, 1, 1;
% 1, 0, 0, 0, 1, 0, 0, 1;
% 0, 0, 0, 0, 0, 0, 1, 1;
% 0, 0, 0, 0, 0, 1, 1, 1;
% 1, 0, 0, 0, 0, 0, 0, 1;
% 1, 1, 0, 0, 0, 0, 1, 1;
% 1, 1, 1, 0, 0, 1, 1, 1 ];
image = [ 1, 0, 1, 0, 1, 0, 1, 0;
0, 1, 0, 1, 0, 1, 0, 1;
1, 0, 1, 0, 1, 0, 1, 0;
0, 1, 0, 1, 0, 1, 0, 1;
1, 1, 1, 1, 1, 1, 1, 1;
0, 0, 0, 0, 0, 0, 0, 0;
1, 1, 1, 1, 1, 1, 1, 1;
0, 0, 0, 0, 0, 0, 0, 0 ];
% Начальный символ приема
sync_symbol = [ 1 ];
% Битовый поток с начальным символом
bit_stream = [ sync_symbol ];
% Генерация битового потока из матрицы профиля pacman
for j = 1:8
bit_stream = [bit_stream image(j,
];end
%bit_stream = (rand(1, 20) > 0.5);
%bit_stream = [1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0];
%bit_stream = [1 1 0 1 0 1 0 1 0];
N = length(bit_stream);
fprintf('bit_stream %d\n', bit_stream);
% Частота семплирования 288 кГц (8 выборок на период)
fs = 288000;
% Несущая частота относительно частоты выборок 36 кГц
f = 0.125;
% Частота доплера 60 Гц (~2.5 м.сек)
fdop = f/600;
% % Частота семплирования 144 кГц (4 выборки на период)
% fs = 144000;
% % Несущая частота относительно частоты выборок 36 кГц
% f = 0.25;
% % Частота доплера 60 Гц (~2.5 м/сек)
% fdop = -f/600;
% Частота среза НЧ фильтра в демодуляторе
fc = 100;
% Частота Найквиста
fn = fs/2;
% Порядок НЧ фильтра демодулятора
lowpfiltord = 50;
% Порядок loop фильтра демодулятора
loopfiltord = 50;
% Амплитуда
A = 1.0;
% Символьная скорость 100 Гц
fsim = 100;
% Количество выборок на период несущей
ns = 1/f;
% Количество выборок необх. для одного символа (360*8 периодов нес. f)
% при символьной скорости fsim
tsymbol = fs*f*ns/fsim;
fprintf('tsymbol %d\n', tsymbol);
% Счетчик выборок в 0
t = 0;
% Количество символов переданных с нулевой фазой (пилот)
nDummySymbols = 4;
% Общее количество выборок
%ttotal = N*tsymbol + nDummySymbols * tsymbol + tsymbol*3;
ttotal = N*tsymbol + nDummySymbols * tsymbol;
fprintf('ttotal %d\n', ttotal);
% Создать пустой BPSK сигнал
BPSK_signal_t = zeros(ttotal,1);
% Пилот сигнал
for i = 1:nDummySymbols
for j = 1:tsymbol
BPSK_signal_t(t+1) = A*sin( 2*pi*(f+fdop)*t );
BPSK_signal_nomod(t+1) = A*sin(2*pi*(f+fdop)*t);
t = t + 1;
end
end
fprintf('t %d\n', t);
% BPSK сигнал
for i = 1:length(bit_stream)
phase = pi * bit_stream(i);
for j = 1:tsymbol
BPSK_signal_t(t+1) = A*sin(2*pi*(f+fdop)*t + phase);
BPSK_signal_mod(t+1) = A*sin(2*pi*(f+fdop)*t + phase);
BPSK_signal_nomod(t+1) = 0;
t = t + 1;
end
end
fprintf('length BPSK_signal %d\n', length(BPSK_signal_t));
fprintf('t %d\n', t);
%% AWGN (аддитивный белый гауссов шум) + полосовая фильтрация
% отношение сигнал/шум в децибелах
SNR = 9;
% мощность сигнала в децибелах
SIGPOWER = 12;
% смесь мод.сигнала с шумом
BPSK_signal_awgn = awgn( BPSK_signal_t, SNR, SIGPOWER );
% Порядок полосового фильтра
bandpfiltord = 50;
% нижняя частота полосового фильтра 35КГц
W1=35e3/(0.5*fs);
% верхняя частота полосового фильтра 37КГц
W2=37e3/(0.5*fs);
% расчет окна Ханна
k=hann(bandpfiltord+1);
% расчет коэффициентов входного полосового фильтра
bp=fir1(bandpfiltord, [W1 W2], k);
% Память состояния фильтра
stateBandpassBpsk = zeros(numel(bp)-1,1);
% Полосовая фильтрация BPSK сигнала + AWGN шум
[BPSK_signal_awgn_filtr,stateBandpassBpsk] = filter(bp,1,BPSK_signal_awgn,stateBandpassBpsk);
% BPSK_signal_awgn_filtr = filter(bp, 1, BPSK_signal_awgn);
% BPSK_signal
%BPSK_signal = BPSK_signal_awgn;
BPSK_signal = BPSK_signal_awgn_filtr;
%% Демодуляция
% полученная картинка в 0
rec_image = zeros(8,8);
% Координаты картинки
x = 1;
y = 1;
% начальный шаг в 0
t = 0;
% Все матрицы в 0
carrier_inph = zeros(ttotal,1);
mixer_inph = zeros(ttotal,1);
demod_thr_inph = zeros(ttotal,1);
%расчет окна Кайзера
kb=kaiser(lowpfiltord+1);
% Коэффициенты НЧ фильтра демодулятора порядок lowpfiltord
b = fir1(lowpfiltord, fc/fn, kb);
% Память состояния для НЧ фильтра демодул.
stateLowpassInph = zeros(numel(
-1,1);% переменные для PLL осциллятора
% cosine output
c = 1;
% delayed cosine by one timestep
c_delay = 0;
% sine output
s = 0;
% delayed sine by one timestep
s_delay = 0;
% Данные из VCO для целей отладки
sine = zeros(ttotal,1);
cosine =zeros(ttotal,1);
vco = zeros(ttotal,1);
%расчет окна Ханна
kp=hann(loopfiltord+1);
% Коэффициенты PLL loop фильтра
bPLL = fir1(loopfiltord, fc/fn, kp);
% Обнулить память
zfPLLa = zeros(numel(bPLL)-1,1);
zfPLLc = zeros(numel(bPLL)-1,1);
zfPLLs = zeros(numel(bPLL)-1,1);
% Обратный отсчет. При нулевом значении отбирается символ.
% Начальный обратный отсчет учитывает задержку
% fir lowpass и пропускает начальный символ
sampleMoment = tsymbol/2+tsymbol;
% Это время необходимое PLL для стабилизации на несущей
nPLLsettle = (nDummySymbols/2) * tsymbol;
% Игнорировать данные, пока не получен символ "1" (синхробит)
ignoreData = 1;
% VCO шаг
VCOgain = 0.005;
% VCO, которое управляет частотой. при v = 0 он точно равен f.
v = 0;
% Основной цикл
for step = 1:ttotal
% Полосовая фильтрация BPSK сигнала + AWGN шум
% [BPSK_signal(step),stateBandpassBpsk] = filter(bp,1,BPSK_signal_awgn(step),stateBandpassBpsk);
% Получить значение BPSK_signal в sample
sample = BPSK_signal(step);
% PLL - управляемый осциллятор с +/- v входом
% w0 частота VCO
w0 = 2*pi*(f+v*VCOgain);
% Cохранить предыдущее состояние
c_delay = c;
s_delay = s;
% Подсчет нового состояния
c = c_delay * cos(w0) - s_delay * sin(w0);
s = s_delay * cos(w0) + c_delay * sin(w0);
% Cохраним все в удобных векторах для графики
sine(step) = s;
cosine(step) = c;
vco(step) = v;
% конец VCO
% Сохраняем sine/cosine. Обратите внимание, что инфазная волна
% теперь является синусоидальной волной, потому что PLL с мультипликацией
% в качестве фазового детектора имеет сдвиг фазы на 90 градусов между
% входом и выходом
carrier_inph(step) = c;
% Это решает, находимся ли мы в режиме синхронизации PLL или режиме замораживания
% Примечание: это должно быть заменено блокировкой PLL в реальном приеме,
% потому что мы не знаем, когда передача фактически
% начинается. Здесь мы просто предполагаем, что она начинается с самого начала
% PLL изменение частоты через v
% Фазовый детектор
pc = sample*c;
ps = sample*-s;
% Фильтрация pc
[vc,zfPLLc] = filter(bPLL,1,pc,zfPLLc);
% Фильтрация ps
[vs,zfPLLs] = filter(bPLL,1,ps,zfPLLs);
[v,zfPLLa] = filter(bPLL,1,vc*vs,zfPLLa);
% шаг на -1
mixer_inph(step) = sample*carrier_inph(step);
[lp_demod_inph(step),stateLowpassInph] = filter(b,1,mixer_inph(step),stateLowpassInph);
demod_thr_inph(step) = lp_demod_inph(step);
if (ignoreData)
% Ждем стартового символа 1
if ( lp_demod_inph(step) > 0.25 )
% Стартовый символ прошел
ignoreData = 0;
end
else
% Ожидание середины символа для выборки демодулятора
sampleMoment = sampleMoment - 1;
if (sampleMoment == 0)
sampleMoment = tsymbol;
if ( demod_thr_inph(step) > 0.1 )
lsb = 1;
else
lsb = 0;
end
fprintf('lsb %d\n', lsb);
% отладка: добавка небольших пиков
lp_demod_inph(step) = lp_demod_inph(step) + 0.1;
% Сохранить результат в картинку 8:8
if (((x<9) && (y<9)) && (ignoreData == 0))
rec_image(y,x) = lsb;
x = x + 1;
if (x > 8)
x = 1;
y = y + 1;
end % if (x > 8)
end % x>9 && y>9
end % sampleMoment == 0
end % ignoreData
% Следующий шаг в выборке
t = t + 1;
end
%% Графика
% % Немодулированный сигнал
% figure;
% subplot(2,1,1);
% plot(BPSK_signal_nomod, 'b-', 'LineWidth', 2);
% hold on;
% grid on;
% title('Немодулированный сигнал (пилот)');
%
% % Модулированный сигнал
% %figure;
% subplot(2,1,2);
% plot(BPSK_signal_mod, 'r-', 'LineWidth', 2);
% hold on;
% grid on;
% title('Модулированный сигнал (данные)');
% Полный BPSK_signal
figure;
subplot(2,1,1);
plot(BPSK_signal_awgn, 'b-', 'LineWidth', 2);
xlabel('Выборки');
hold on;
grid on;
title('Полный BPSK сигнал + шум AWGN');
subplot(2,1,2);
plot(BPSK_signal, 'r-', 'LineWidth', 2);
xlabel('Выборки');
hold on;
grid on;
title('Полный BPSK сигнал + шум AWGN + Фильтр');
% Демодулированный сигнал
figure;
subplot(2,1,1);
plot(lp_demod_inph,'LineWidth',2);
xlabel('Выборки');
ylabel('Амплитуда');
title('Демодулированный сигнал / НЧ фильтрация');
grid on;
% Картинка
subplot(2,1,2);
imshow(rec_image);
title('Битовая картинка 8:8 квадратиков 1-белый, 0-черный кв. ' );
% Напряжение VCO
% figure
% plot(vco);
% title('Напряжение VCO');
%
% % Фильтр демодулятора
% figure
% % АЧХ НЧ фильтра
% [h,f]=freqz(b,1);
% % Параметры являются соответственно частотой и амплитудой
% plot(f*fs/(2*pi),20*log10(abs(h)))
% xlabel('частота/Hz');ylabel('амплитуда/dB');title('АЧХ НЧ фильтра демодулятора');
%
% % Фильтр полосовой на 36 кГц
% figure
% [h,fp]=freqz(bp,1);%amplitude-frequency characteristic graph
% plot(fp*fs/(2*pi),20*log10(abs(h)))%parameters are respectively frequency and amplitude
% xlabel('частота/Hz');ylabel('амплитуда/dB');title('АЧХ полосового фильтра');
Русская версия Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)