Всем привет.
Нужно создать сигнал который уже находится на выходе квадратурного демодулятора, то есть его огибающую. Способ модуляции 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 это не понятно. Да и ко всему она же каждый раз меняется.
Может кто знает в чем моя ошибка?
Заранее спасибо.

Вы построили периодограмму одной реализации случайного процесса. Она, естественно, случайна. Чтобы получить оценку спектра периодограммы надо усреднять.
Пока усредняете, подумайте, какой вообще должен получиться спектр и какому отсчету у вас будет соответствовать частота 2.5 кГц.

пасиб попробуем. есть в матлабе функция "periodogram" ... что то в этом роде... сейчас ею займусь
...............
незнаю или я правильно сделал, но после строчки
msg_filt = filter(coefs,1,s_psk_300);сглаживание углов
я добавил еще две
carrier = exp(1i*15.7e3*t);
y=msg_filt.*carrier;

как бы сдвинул спектр. незнаю только почему (15.7 это 2п*2.5е3) нужная частота не получилась сразу.
может если разберусь с "periodogram" тогда и получу ответ
Често говоря, в этой теме не силен... но мне кажется что случайность информации не должна влиять на частоту модуляции...ведь частота модуляции она детерминирована

Рекомендую пока убрать фильтрацию и сформировать сигнал с прямоугольными символами. Для такого сигнала посчитайте АКФ, а из нее СПМ, будет ясно, почему нету максимума в 2.5 кГц.

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

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

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



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



По большому счету да. Но при моих условиях они совпадают (почти). для большей наглядности я взял 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кГц тоже переместился на соответствующую величину.
Вот такие пироги

Но АКФ действительная. Если есть сложности с комплексным сигналом, возмите ФМ2 с огибающей +-1. Результат будет тот же. Решите эту задачу на бумаге прежде, чем что-то моделировать в матлабе. Если совсем уж туго пойдет, подсмотрите в книжке. У Прокиса, например.

Теперь про ваш код.
В предыдущем варианте результат больше был похож на правду.
Не знаю, как работает функция periodogramm. Рекомендую писать в лоб: разбить реализацию на окна (функция buffer), взять от каждого ДПФ, и усреднить квадрать модуля ДПФ по всем окнам. Выглядит это примерно так: spectrum=mean(abs(fft(buffer(signal,window_size))).^2, 2);

Сначала получите спектр на бумаге, потом получите такой же в матлабе. Только после этого имеет смысл обсуждать положение максимума.

Только в нуле. Оно и понятно - в нуле мы получаем power.
Более того... для того чтоб доказать что АКФ для недействительного сигнала недействилелна, нужно просто убедится что
conj(Rxx[k]) not= Rxx[k]. Взяв за основу определение АКФ для комплексных сигналов, в этом убедиться не сложно.
Лично я пошел дальше и для себя доказал что conj(Rxx[-k]) = Rxx[k]. Уверенн что это было доказано уже миллион раз но мне
быстрее было доказать самому нежели искать в итернете или книжках.
В общем последняя строчка (которую я написл выше) мне совсем не понравилась... особенно для ДПФ...

В плане оценки частот - да. Но не в плане АКФ и ДПФ от АКФ нет

Так в том то и дело что в свое время перерешал таких задач много и все на бумаге. А вот теперь пошла практика. Все мы прекрасно знаем разницу меду обеими (:

Вы имеете в виду мое умножение на экспоненту?
Честно говоря я тут тоже сомневаюсь. Ведь умножение на эксп. это всеравно что умножить на син. и кос. (по Ойлеру)
Тоесть я поменял "природу" самого сигнала...верно...
может вместо умножения на експ. нужно сделать что то со временем. Типа scaling ?

То что вы предлагаете это называется Bartlett method
http://en.wikipedia.org/wiki/Bartlett%27s_method
В виду того что я взял много выборок мне уже не нужно оценивать PSD. Она у меня уже явно детерминированная (судя по стабильности рузультатов). Возьмите количество символов 75 и тогда будут видны (совсем небольшие) изменения.

Везде.


Во всем такой же результат будет.


Ну так какой спектр должен быть в теории?


Нет. Ваш спектр не похож на тот, который должен быть.



Да. Вы перенесли сигнал на несущую.


Ничего такого не нужно.


Ну да. И он с ходу дает то, что должно быть.

Если не тяжело - докажите или хотя бы дайте ссылку.
Я ведь обосновал почему я думаю что она коплексная а вы просто написали одно слово.


Зависит от "Докажите или дайте ссилку"...


Ассиметричный если сигнал действительный (да еще и комплексный).
А есле сигнал еще и комлексный то какой спектр должен быть?


А какой должен быть?
Я всего навсего изменил количество символов


Мне кажется что тут надо подумать.Ведь я по большому счету нигде не указал что мой
сигнал "живет" на частоте 2.5кГц. Единственное что (судя из моего кода и результатов отклика Матлаба) я указал четко так это частоту дискретизации.
По этому с вашего позволения я над этим подумаю


Зачем мне оценка PSD если у меня есть настоящий?

Корреляция между отсчетами сигнала в любые два момента времени будет равна либо среднему квадрату амплитуды, либо нулю, если символы некоррелированы.
Книжку я указал: Прокис "Цифровая связь".


У него аналитическое выражение есть.


Посчитайте уже, наконец, или прочитайте.


Еще раз. У вас - только оценка СПМ, настоящая будет только на бумаге.

Предлагаю решать по одному вопросу. Первым делом разберемся с АКФ комплексного сигнала.

Я просто нехочу повторятся. Поэтому если не тяжело дайте ссылку на книгу и только подскажите в какой главе
главе говориться про автокореляцию комплексного сигнала. Я всего навсего хочу убедится что она действительна. Хотя ссылки на книгу уже будет достаточно.

Более того... сдесь http://en.wikipedia.org/wiki/Autocorrelation под заголовком Properties Вы обнаружите что она комплексная.
Как я уже говорил поначалу я это доказал сам исходя из определения АКФ

Не вижу проблемы. Действительное число - частный случай комплексного.

Вот именно... частный. А мы говорим об общем.

В вашем частном случае АКФ действительна.

А разве сигнал у меня действительный?

Нет. Его АКФ - да.

Так я же ссылку Вам показывал.... Там же написано было что АКФ комплексного сигнала комплексная

Ну ладно, комплексная. С нулевой мнимой частью.

A fundamental property of the autocorrelation is symmetry, R(i) = R(-i), which is easy to prove from the definition. In the continuous case,

the autocorrelation is an even function

R_f(-tau) = R_f(au), when f is a real function

and the autocorrelation is a Hermitian function

R_f(-tau) = R_f^*(tau), when f is a complex function.

Это выдержка из Вики (ссылки которую я вам давал)
1)Где тут написано что мнимая часть априори равна нулю?
2)Зачем применять conj к выражению у которого мнимая часть априори равна нулю (тоесть к действительному выражению)?

Как это все противоречит написанному мной?
Учите математику, читайте книжку. Больше пока ничем помочь не могу.

Как? Да на прямую. По поводу "мнимой нулевой" части
Я и так знаю что мне нужно делать. Судя по Вами вышесказанному математику я знаю не хуже Вас.
Читать книжки это полезно... но лучше "иногда" соображать самому.
Мы немного начинаем уходить от темы.
Давайте другим тоже дадим высказаться (по моему вопросу который я задал в самом начале)