clearvars %------------------------------------------------------------------------ % Задаём параметры сигналов и поток битов bits = [1 0 1 1 1 0 0 1 0 0 1 0 1 1 1 0 0 1 0 1 randi(2, [1,9980])-1]; % поток битов N = length(bits); % число битов R = 5*10^6; % количество передаваемых бит в единицу времени (бит/с) T = 1 / R; % длительность одного бита dev_f = R / 2; % девиация частоты при передаче битов для условия их ортогональности при некогер. приёме. % Интервал частот между символами равен R, тогда девиация частоты - R/2. % Общая эффективная ширина спектра сигнала - 2R. Pn = 0.3; % Мощность БГШ в Вт на нагрузке 1 Ом. f0 = 50 * 10^6; % частота несущей, Гц Fs = 10 * f0; % частота дискретизации Ts = 1/Fs; % Период дискретизации FsFd = round(Fs * T); % число временных отсчётов в одном бите Signal_duration = FsFd * N; % Число временных отсчётов в исходном сигнале, длительность ТОЛЬКО СИГНАЛА t = (0:3*Signal_duration/2-1)*Ts; % временной вектор всего потока данных Duration = length(t); % число временных отсчётов в реализации, длительность ВСЕЙ РЕАЛИЗАЦИИ f = Fs/2 * linspace(0, 1, Duration/2+1); % вектор отсчётов частоты: от 0 до 1/2 частоты дискретизации % Количество отсчётов БПФ совпадает с кол-ом отсчётов в реализации. % Мы берём половину отсчётов БПФ. disp('Длительность одного бита') disp(T) load('Coeff_A_file.mat') % загружаем коэффициенты фильтра. Коэффициенты в переменной Coeff_A. F_Delay = 89; % Групповая задержка появления сигнала на выходе фильтра в отсчётах времени % Эта переменная нужна для учёта сдвига во времени выходного сигнала % фильтра. Она используется при вычислении мощности сигнала и в приёмнике % для выделения интервалов времени интегрирования. Av = 1; % Амплитуда видеоимпульсов передатчика, В Ac = 1; % Амплитуда колебания несущей, В %------------------------------------------------------------------------ % ММС сигнал % Формирование через формулу косинуса суммы двух углов % Формируем видеосигнал из потока битов Video_Pulse_Tr = zeros(1, Duration); for i = 1 : N if bits(i) == 0 for j = (i-1)*FsFd+1 : i*FsFd Video_Pulse_Tr(j) = -Av; end else for j = (i-1)*FsFd+1 : i*FsFd Video_Pulse_Tr(j) = Av; end end end % Строим график видеосигнала % figure % plot(t, Video_Pulse_Tr) % grid on % Формируем модулирующее изменение полной фазы несущей % Формирование проводим только для тех отсчётов времени, когда есть поток % битов Phi = 2* pi * dev_f * Ts* cumsum(Video_Pulse_Tr(1:Signal_duration)); % Строим график этого сигнала % figure % plot(t(1:5*FsFd), 180/pi*Phi(1:5*FsFd)) % grid on % Формируем компоненты для модуляции несущей % Они тоже формируются только на интервале времени потока битов I = cos(Phi); Q = sin(Phi); % Строим графики % figure % plot(t(1:5*FsFd+1), I(1:5*FsFd+1), t(1:5*FsFd+1), Q(1:5*FsFd+1)) % grid on % Формируем синфазную и квадратурную компоненты несущей % На интервале потока битов. I_carr = Ac * cos(2*pi*f0*t(1:Signal_duration)); Q_carr = Ac * sin(2*pi*f0*t(1:Signal_duration)); % Формируем ЧМнФ сигнал % На интервале потока битов S = I.*I_carr - Q.*Q_carr; % Дополняем сигнал нулями для учёта фильтра в дальнейшем % Дополняем до временного интервала всей реализации S = [S zeros(1,Duration-Signal_duration)]; % Вычисляем мощность ММС сигнала до формирующего фильтра, нагрузка 1 Ом % Dens_Before = var( S(1:Signal_duration) ); % Пропускаем ММС-сигнал через фильтр S1 = filter(Coeff_A, 1, S); % Вычисляем мощность после фильтра Dens_After = var( S1(1+F_Delay:Signal_duration+F_Delay) ); % Строим осциллограмму ММС до фильтра % figure % plot(t(1:5*FsFd+1), S(1:5*FsFd+1)) % grid on % Строим осциллограмму ММС после фильтра % figure % plot(t, S1) % grid on % grid minor %------------------------------------------------------------------------- % % Вычисляем СПМ получившегося сигнала до фильтра % S_f = fft(S)/Duration; % вычисляем БПФ % psd_S_f = S_f.*conj(S_f); % вычисляем квадрат модуля БПФ, т.е. СПМ % psd_S_f1 = psd_S_f(1:Duration/2+1); % выделяем из всей СПМ отсчёты от 0 до 1/2 частоты дискретизации % psd_S_f1(2:Duration/2) = 2 * psd_S_f1(2:Duration/2); % умножаем выделенные отсчёты на 2, % % чтобы учесть энергию той части спектра, которую мы убрали из % % рассмотрения, нулевой и последний отсчёты на 2 не умножаем. % % Вычисляем СПМ получившегося сигнала после фильтра % S1_f = fft(S1)/Duration; % вычисляем БПФ % psd_S1_f = S1_f.*conj(S1_f); % вычисляем квадрат модуля БПФ, т.е. СПМ % psd_S1_f1 = psd_S1_f(1:Duration/2+1); % выделяем из всей СПМ отсчёты от 0 до 1/2 частоты дискретизации % psd_S1_f1(2:Duration/2) = 2 * psd_S1_f1(2:Duration/2); % умножаем выделенные отсчёты на 2, % % чтобы учесть энергию той части спектра, которую мы убрали из % % рассмотрения, нулевой и последний отсчёты на 2 не умножаем. % Строим графики СПМ % figure % plot(f, 10*log10(psd_S_f1*10^3)) % xlabel('f, Гц') % ylabel('PSD, дБ (мВт / Гц)') % title('СПМ ММС-сигнала до и после фильтра') % % legend('До фильтра','После фильтра') % set(gca, 'Ylim', [-100, 30], 'Xlim', [5*10^6, 95*10^6]) % grid on % grid minor %-------------------------------------------------------------------- % Создаём реализацию белого шума WNoise = wgn(1, Duration, Pn, 'linear', 'real'); % шум, Pn (Вт) - мощность белого шума (дисперсия), нагрузка 1 Ом % BandNoise = filter(Coeff_A, 1, WNoise); % фильтруем БГШ Sigma = var(WNoise); % вычисляем мощность получившегося шума SNR = Dens_After / Sigma; % вычисляем сигнал/шум по мощности, нагрузка 1 Ом Process_in = S1 + WNoise; % складываем ММС-сигнал и белый ГШ % % Вычисляем СПМ получившейся аддитивной смеси % Process_in_f = fft(Process_in)/Duration; % вычисляем БПФ % psd_Process_in_f = Process_in_f.*conj(Process_in_f); % вычисляем квадрат модуля БПФ, т.е. СПМ % psd_Process_in_f1 = psd_Process_in_f(1:Duration/2+1); % выделяем из всей СПМ отсчёты от 0 до 1/2 частоты дискретизации % psd_Process_in_f1(2:Duration/2) = 2 * psd_Process_in_f1(2:Duration/2); % умножаем выделенные отсчёты на 2, % % чтобы учесть энергию той части спектра, которую мы убрали из % % рассмотрения, нулевой и последний отсчёты на 2 не умножаем. % % % Строим график СПМ % figure % plot(f, 10*log10(psd_Process_in_f1*10^3)) % xlabel('f, Гц') % ylabel('PSD, дБ (мВт / Гц)') % title('СПМ суммы ММС-сигнала и полосового ГШ, вход приёмника') % set(gca, 'Ylim', [-100, 30], 'Xlim', [5*10^6, 95*10^6]) % grid on % grid minor %------------------------------------------------------------------------ % Пропускаем аддитивную смесь сигнала и шума через ограничитель % Process_in_1 = zeros(1, Duration); % for i = 1 : Duration % if Process_in(i) > 0 % Process_in_1(i) = 1; % else % Process_in_1(i) = 0; % end % end % Фильтруем основную гармонику несущего колебания на выходе % ограничителя % Process_out = filter(Coeff_A, 1, Process_in_1); % % Вычисляем СПМ процесса на выходе ограничителя % Process_out_f = fft(Process_out)/Duration; % вычисляем БПФ % psd_Process_out_f = Process_out_f.*conj(Process_out_f); % вычисляем квадрат модуля БПФ, т.е. СПМ % psd_Process_out_f1 = psd_Process_out_f(1:Duration/2+1); % выделяем из всей СПМ отсчёты от 0 до 1/2 частоты дискретизации % psd_Process_out_f1(2:Duration/2) = 2 * psd_Process_out_f1(2:Duration/2); % умножаем выделенные отсчёты на 2, % % чтобы учесть энергию той части спектра, которую мы убрали из % % рассмотрения, нулевой и последний отсчёты на 2 не умножаем. % Строим график СПМ на выходе % figure % plot(f, 10*log10(psd_Process_out_f1*10^3)) % xlabel('f, Гц') % ylabel('PSD, дБ (мВт / Гц)') % %title('СПМ суммы ММС-сигнала и полосового ГШ, вход приёмника') % set(gca, 'Ylim', [-100, 30], 'Xlim', [5*10^6, 95*10^6]) % grid on % grid minor %------------------------------------------------------------------------- % Реализуем некогерентный приём ММС-сигнала % Начальные фазы гетеродинов в обои каналах выбраны произвольно, т.к. % некогерентный приём. %Канал выделения символа 0 %------------------------- % Формируем квадратуры гетеродина частоты символа 0 LVCO_0_cos = cos(2*pi*(f0-dev_f)*t+0.12*pi); % косинус гетеродина LVCO_0_sin = sin(2*pi*(f0-dev_f)*t+0.12*pi); % синус гетеродина % Умножаем каждую квадратуру с пришедшим сигналом Receive10 = Process_in .* LVCO_0_cos; % косинус Receive20 = Process_in .* LVCO_0_sin; % синус % Интегрируем произведения и возводим результат в квадрат на интервале передачи одного бита Receive30 = zeros(1, Duration); % массив результатов интегрирования и квадрата для косинуса Receive40 = zeros(1, Duration); % массив результатов интегрирования и квадрата для синуса for i = 1 : N Bit_interval = Receive10( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ); % выделяем один символ в косинусе Rv1 = Ts * cumsum(Bit_interval); % интегрируем Rv1 = Rv1 .* Rv1; % возводим в квадрат Receive30( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ) = Rv1; % записываем в массив результатов для косинуса Bit_interval = Receive20( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ); % выделяем один символ в синусе Rv1 = Ts * cumsum(Bit_interval); % интегрируем Rv1 = Rv1 .* Rv1; % возводим в квадрат Receive40( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ) = Rv1; % записываем в массив результатов для синуса end % Складываем их Receive50 = Receive30 + Receive40; % Извлекаем квадратный корень из суммы Receive60 = sqrt(Receive50); %Канал выделения символа 1 %------------------------- % Формируем квадратуры гетеродина частоты символа 1 LVCO_1_cos = cos(2*pi*(f0+dev_f)*t+0.61*pi); % косинус гетеродина LVCO_1_sin = sin(2*pi*(f0+dev_f)*t+0.61*pi); % синус гетеродина % Умножаем каждую квадратуру с пришедшим сигналом Receive11 = Process_in .* LVCO_1_cos; % косинус Receive21 = Process_in .* LVCO_1_sin; % синус % Интегрируем произведения на интервале передачи одного бита Receive31 = zeros(1, Duration); % массив результатов для косинуса Receive41 = zeros(1, Duration); % массив результатов для синуса for i = 1 : N Bit_interval = Receive11( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ); % выделяем один символ в косинусе Rv1 = Ts * cumsum(Bit_interval); % интегрируем Rv1 = Rv1 .* Rv1; % возводим в квадрат Receive31( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ) = Rv1; % записываем в массив результатов для косинуса Bit_interval = Receive21( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ); % выделяем один символ в синусе Rv1 = Ts * cumsum(Bit_interval); % интегрируем Rv1 = Rv1 .* Rv1; % возводим в квадрат Receive41( (i-1)*FsFd+1+F_Delay : i*FsFd+F_Delay ) = Rv1; % записываем в массив результатов для синуса end % Складываем их Receive51 = Receive31 + Receive41; % Извлекаем квадратный корень из суммы Receive61 = sqrt(Receive51); %Общий канал %----------- % Вычитаем из канала 0 канал 1 Receive7 = Receive60 - Receive61; % Из получившегося процесса выбираем последние отсчёты в каждом % символе, т.е. 40-ой отсчёт, 80-й, 120-й и т.д. Solve_data = zeros(1, N); % массив этих отсчётов for i = 1 : N Solve_data(i) = Receive7(i*FsFd+F_Delay); end % Выносим решение о том, какой бит был передан. Если i-ый элемент % массива >0, был передан бит 0, если <0, был передан бит 1. Received_bits = zeros(1, N); for i = 1 : N if Solve_data(i) > 0 Received_bits(i) = 0; else Received_bits(i) = 1; end end % Смотрим, сколько ошибок произошло Error = 0; for i = 1 : N if bits(i) ~= Received_bits(i) Error = Error + 1; end end disp('Число ошибок') disp(Error) disp('Отношение сигнал / шум на входе приёмника, раз') disp(SNR)