Цитата(ataradov @ Apr 29 2018, 18:35)
Чем обычная fft() не устраивает? Нужно только результат в dB перевести, но это тривиально.
В моем наборе функций есть fftplot(). Она довольно топорная, но посмотреть можно.
fft() почему-то рисует немного не то. Спектр смещен вправо по оси частот. Не врубаюсь, почему?
CODE
clear all;
close all;
%% Исходные установки
% Битовая посл.
%bit_stream = (rand(1, 20) > 0.5);
bit_stream = [0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
%bit_stream = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% длина бит. посл.
N = length(bit_stream);
% Несущая частота Гц
f = 36000;
% Кол. выб. на период несущей
nsp = 4;
% Частота семплирования Гц
fs = f*nsp;
% Частота доплера Гц
fdop = 0;
% время выборки
Ts = 1/fs;
fprintf('Ts %d\n', Ts);
% Амплитуда
A = 0.1;
% колич. периодов несущей для одного симв.
sone = 360;
% Символьная частота Гц
fsim = f/sone;
fprintf('fsim %d\n', fsim);
% Количество выборок необх. для одного символа
% при символьной скорости fsim
tsymbol = f*nsp/fsim;
fprintf('tsymbol %d\n', tsymbol);
% несущая частота для PLL Гц
fref = f;
% Доплер для PLL
ferr = fdop;
fprintf('fref %d\n', fref);
ttotal = N*tsymbol;
% Частота Найквиста
fn = fs/2;
% несущая и доплер в % по отн. к fs
fb = f/fs;
fbdop = fdop/fs;
%% Входной сигнал
t = 0;
for i = 1:N
% phase = pi * bit_stream(i);
phase = 0;
for j = 1:tsymbol
% входной сигнал
BPSK_signal(t+1) = cos(2*pi*(fb+fbdop)*t + phase);
t = t + 1;
end
end
%% Графика
figure
psd(BPSK_signal(100:end),2^14,fs/1e3)
axis([35 37 -120 40]),xlabel('КHz'),title('Спектр сигнала')
figure
Nb = length(BPSK_signal);
X_mags = abs(fft(BPSK_signal));
bin_vals = [10 : Nb-1];
fax_Hz = bin_vals*fs/Nb;
Nb_2 = ceil(Nb/2);
plot(fax_Hz(1:Nb_2), 10*log10(X_mags(1:Nb_2)))
axis([35.5e3 36.5e3 -120 40]),xlabel('КHz'),title('Спектр сигнала')
Нажмите для просмотра прикрепленного файла Нажмите для просмотра прикрепленного файлаА, разобрался.
Спасибо.
Примерно так:
Код
figure
Nb = length(BPSK_signal);
X_mags = abs(fft(BPSK_signal));
bin_vals = [1 : Nb];
fax_Hz = bin_vals*fs/Nb;
plot(fax_Hz(1:Nb), 10*log10(X_mags(1:Nb)))
axis([35.5e3 36.5e3 -120 40]),xlabel('Hz'),title('Спектр сигнала')