Создать выборки огибающей QPSK сигнала в Матлабе, на выходе квадратурного демодулятора Опции

[reginil_y]

Всем привет.
Нужно создать сигнал который уже находится на выходе квадратурного демодулятора, то есть его огибающую. Способ модуляции QPSK.
Количество выборок равно 300, 4 выборки на символ. Получается всего 75 символов. Символьная скорость (Fd) равна 2.5кHz. Частота дикретизации равна 10кHz (FS).
привожу код:

%Постройка огибающей
Nsym = 75;
Fd = 2.5e3;
Fs = 10e3;
coefs = rcosine(Fd,Fs);%подготовка коэфициентов для формирующего фильтра
msg = [randint(1, Nsym, 3)]; % подготовка символов
t = (0 : 1/Fs : 1/Fd*Nsym-1/Fs); % дискретное время
s_qpsk = pskmod(msg, 4,pi/4); % собственно QPSK сигнал
s_psk_300 = s_qpsk(floor(Fd*t)+1); % постройка вектора длинной в 300. (4 выборки на символ)
msg_filt = filter(coefs,1,s_psk_300);сглаживание углов

а теперь строю периодограмму и смотрю на какой частоте (fec) получается максимум (желательно было бы на 2.5 кHz)

%Собственно построение периодограммы
resf=2^14;
fftr = abs(fft(msg_filt,resf)/sqrt(resf)).^2;
figure
plot(0:1/resf:1-1/resf,fftr)
title('FFT of msg-filt')

[max1 ind1] = max(fftr);
ind1=ind1-1;
if ind1<resf/2
fec = (ind1/(resf/Fs))
else
fec = (-(resf-ind1)/(resf/Fs))
end

То что спектр не симетричен это понятно - ДПФ комплекного сигнала.
А вот почему оцениваемая частота не равна Fd это не понятно. Да и ко всему она же каждый раз меняется.
Может кто знает в чем моя ошибка?
Заранее спасибо.

[KalashKS]

АКФ делать не стал по двум причинам. 1-ая: тажело будет глядя на функцию во временной области понять на глаз что с ней произойдет в частотной области и почему пик не выскочит в 2.5 кило герц. ну а вторая причина это то что СПМ это и есть периодограмма. Если хотите могу разместить код.

Понять будет легко. Формы что АКФ, что спектра простые. Ваша проблема в том, что не зная теории, вы не можете отличить правильный результат от неправильного.
И да. Периодограмма - не то же самое, что СПМ, а только ее оценка. Более того, оценка смещенная.
Код - размещайте.

[reginil_y]

Понять будет легко. Формы что АКФ, что спектра простые.

Я сомневаюсь, так как огибающая комплексная

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

может быть да... мне немного нехватает теории... это тоже одна из причин по которой я разместил вопрос

И да. Периодограмма - не то же самое, что СПМ, а только ее оценка. Более того, оценка смещенная.

По большому счету да. Но при моих условиях они совпадают (почти). для большей наглядности я взял Nsym=2000. Тоесть количество выборок равно 8000. Ниже я размещу код который убедит вас в том что несмотря на случайность процесса результаты "оценки" частот абсолютно стабильны! Он так же вас убедит в том что формирующий фильтр ничего не меняет.

Правда осталось пару детерминированных вопросов.
Вот код:

close all
clear all
clc

resf=2^15;
Nsym = 2000;
Fd = 2.5e3; % modulation frequency
Fs = 10e3; % sampling frequency
coefs = rcosine(Fd,Fs,'fir/sqrt');
msg = [randint(1, Nsym, 3)]; % creating the massage from 75 symbols
t = (0 : 1/Fs : 1/Fd*Nsym-1/Fs); % descrete time
s_qpsk = pskmod(msg, 4,pi/4); % the modulation signal
s_psk_300 = s_qpsk(floor(Fd*t)+1); % lifting up the descrite time

[Pxx_without_f,w_without_f] = periodogram(s_psk_300,[],resf);

figure(1)
plot(w_without_f,Pxx_without_f)
title('Periodogram (Matlab functioin) before filter ')
[max1m ind1m] = max(Pxx_without_f);
ind1m=ind1m-1;
if ind1m<w_without_f/2
fecm = (ind1m/(resf/Fs));
else
fecm = (-(resf-ind1m)/(resf/Fs));
end

figure(2)
plot(real(s_psk_300),'.-r')
hold on
plot(imag(s_psk_300),'.-b')
title('The Signal before Filter')
axis([-2 100 -2 2])

msg_filt = filter(coefs,1,s_psk_300);

[Pxx_with_f,w_with_f] = periodogram(msg_filt,[],resf);
figure(3)
plot(w_with_f,Pxx_with_f)
title('Periodogram (Matlab functioin) after filter ')
[max1a ind1a] = max(Pxx_with_f);
ind1a=ind1a-1;
if ind1a<w_with_f/2
feca = (ind1a/(resf/Fs));
else
feca = (-(resf-ind1a)/(resf/Fs));
end



carrier = exp(-2i*pi*2.5e3*t);
y=msg_filt.*carrier;

[Pxx_after_shift,w_after_shift] = periodogram(y,[],resf);
figure(4)
plot(w_after_shift,Pxx_after_shift)
title('Periodogram (Matlab functioin) after shift ')
[max1s ind1s] = max(Pxx_after_shift);
ind1s=ind1s-1;
if ind1s<w_after_shift/2
fecs = (ind1s/(resf/Fs));
else
fecs = (-(resf-ind1s)/(resf/Fs));
end


figure(5)
plot(real(msg_filt),'.-r')
hold on
plot(imag(msg_filt),'.-b')
title('Output of the shapinfg filter')
axis([-2 100 -5 5])

figure(4)
plot(real(y),'.-r')
hold on
plot(imag(y),'.-b')
title('y After shift to 2.5e3Hz')
axis([-2 100 -5 5])


%Creating PSD for original signal
fftr = abs(fft(s_psk_300,resf)/sqrt(resf)).^2;
figure
plot(0:1/resf:1-1/resf,fftr)
%axis([-0.1 1.1 -0.1 14])
title('FFT of msg-filt')

[max1 ind1] = max(fftr);
ind1=ind1-1;
if ind1<resf/2
fec = (ind1/(resf/Fs));
else
fec = (-(resf-ind1)/(resf/Fs));
end

%Creating PSD for the signal after filtering
fftr1 = abs(fft(msg_filt,resf)/sqrt(resf)).^2;
[max11 ind11] = max(fftr1);
ind11=ind11-1;
if ind11<resf/2
fec1 = (ind11/(resf/Fs));
else
fec1 = (-(resf-ind11)/(resf/Fs));
end

%Creating PSD for the signal after filtering
fftr2 = abs(fft(y,resf)/sqrt(resf)).^2;
[max12 ind12] = max(fftr2);
ind12=ind12-1;
if ind12<resf/2
fec2 = (ind12/(resf/Fs));
else
fec2 = (-(resf-ind12)/(resf/Fs));
end


disp([' Estimation frec. before filt. = ' num2str(fecm) ' Estimation frec. after filt. = ' num2str(feca) ' Estimation frec. after shift = ' num2str(fecs)])
disp([' Estimation frec. before filt.(my func.) = ' num2str(fec)])
disp([' Estimation frec. after filt.(my func.) = ' num2str(fec1)])
disp([' Estimation frec. after shit(my func.) = ' num2str(fec2)])

вот (стабильный) отклик матлаба на него

Estimation frec. before filt. = -10000 Estimation frec. after filt. = -10000 Estimation frec. after shift = -2500
Estimation frec. before filt.(my func.) = 0
Estimation frec. after filt.(my func.) = 0
Estimation frec. after shit(my func.) = -2500

теперь вместо carrier = exp(-2i*pi*2.5e3*t); делаем carrier = exp(2i*pi*2.5e3*t);

получаем немного другой но тоже стабильный отклик

Estimation frec. before filt. = -10000 Estimation frec. after filt. = -10000 Estimation frec. after shift = -7500
Estimation frec. before filt.(my func.) = 0
Estimation frec. after filt.(my func.) = 0
Estimation frec. after shit(my func.) = 2500

То есть я каждый раз вычисляю при какой частоте получу макс. значение периодограммы и abs(fft(.)/N)^2... то что я назвал my func.

Прежде всего видно что формирующий фильтр ни на что не повлиял как в периодограмме так и в my func. (уже хорошо)

первый вопрос (главный) который возникает сам собой - почему макс. значение периодограммы в -10 кГц ..... (и вообще - что такое -10 кГц) а в my func. в нуле (может это одно и тоже?)

Ну а второй вопрос по поводу сдвижки спектра
В случае с периодограммой: при максимуме в точке -10кГц и при сдвижке спектра на -2.5кГц получаю максимум в -2500Гц. то есть вроде как действительно -10кГц это как 0Гц. А при сдвижке спектра на 2.5кГц получаю максимум в -7500Гц. что более логично чем в первом случае.
При работе с my func. все как то более понятно (вопрос правиль но ли?)
в обеих случаях максимум имел место в нуле. полсле сдвижки на -2.5кГц он переместился именно на эту величину а при сдвижке на 2.5кГц тоже переместился на соответствующую величину.
Вот такие пироги