День 3 Летней школы Julia
Часть 1: генерация сигналов
Рассмотрим основные техники генерации периодических, импульсных и шумовых сигналов для тестирования алгоритмов цифровой обработки сигналов (ЦОС) и радиотехнических систем.
Используемые библиотеки:
#Pkg.add("DSP")
#Pkg.add("ChirpSignal")
#Pkg.add("Waveforms")
using DSP, ChirpSignal, Waveforms
Генерация периодических сигналов
Периодический сигнал - это повторяющаяся последовательность. Создадим один период значений дискретного сигнала и визуализируем во времени:
one_period = [1; 4; 6; 8; 7; 5; 0; 0; 0; 0];
plot(one_period, line=:stem, linewidth=2, marker=:circle)
Повторим период пять раз функцией repeat. Визуализируем во времени ступенчатой линией:
five_period = repeat(one_period, 5);
plot(five_period, linewidth=2, marker=:circle, line=:steppost)
Попробуем сгенерировать периодический сигнал при помощи математических функций, таких, к примеру, какsin. Для того, чтобы сгенерировать синусоиду с требуемой частотой повторения, стоит воспользоваться вектором времени - последовательностью равномерно нарастающих отсчётов времени с указанным периодом дискретизации. Зададим вектор времени для сигналов с частотой дискретизации 1000 Гц:
fs = 1000; # частота дискретизации
dt = 1/fs; # шаг временной сетки
stoptime = 0.5;
t = 0:dt:stoptime-dt
Теперь создадим периодический синусоидальный сигнал с основной частотой 10 Гц (форма синусоиды будет повторяться каждые 0.1 сек). Переменная sine_wave - это вектор из 500 отсчётов.
sine_wave = sin.(2*pi*10*t);
plot(t, sine_wave)
Мы можем одной строчкой кода сгенерировать несколько синусоид с разными амплитудами, частотами и начальными фазами. Для этого зададим эти параметры векторами. Результирующая переменная three_sines уже будет матрицей размером 500х3. Мы можем передать её функции plot для одновременного отображения трёх графиков на одних осях:
three_sines = [1.4 0.6 1] .* sin.(2*pi*t.*[10 15 30] .+ [pi/4 0 pi/6]);
plot(t, three_sines)
Если нам интересно наблюдать форму суммы трёх синусоид, то мы можем сложить столбцы матрицы при помощи функции sum:
sum_sines = sum(three_sines, dims = 2);
plot(t, sum_sines, title = "Сумма трёх синусоид")
Дополнительные функции для генерации сигналов
Мы можем генерировать одни периодические сигналы на основе других. Например, получать форму типичных пилообразных и треугольных сигналов. Сгенерируем меандр из синусоиды:
meandr = sine_wave .> 0;
plot(t, meandr)
Процент "заполнения" периода сигнала, или его скважность, можно легко контролировать пороговым уровнем:
square_wave = sine_wave .> -0.5;
plot(t, square_wave)
Пилообразный сигнал, изменяющийся в пределах от 0 до +1 можно получить функцией mod, передав ей вектор времени:
sawtooth = mod.(t*10,1);
plot(t, sawtooth)
Библиотека Waveforms.jl содержит функции для генерации различных периодических сигналов, но в качестве "задающего" сигнала, определяющего основную частоту, используется равномерно нарастающий пилообразный сигнал, изменяющийся в пределах от 0 до 2pi. См. документацию.
Случайные величины и шумовой сигнал
В Engee доступны функции для генерации случайных чисел с различными распределениями. Самые популярные - это равномерное и нормальное. Рассмотрим применение встроенных функций rand и randn:
rand()
rand(3)
rand(1:100)
rand(1:100, 6)
Воспользуемся функцией rand для создания шумового сигнала с равномерным распределением:
noise = rand(length(t));
noise = noise .- 0.5;
plot(t, noise)
Убедимся, что шум с равномерным распределением, построив гистограмму вектора:
histogram(noise; nbins = 20)
Шум с Гауссовским распределением:
noise_norm = randn(length(t));
plot(t, noise_norm)
И его гистограмма:
histogram(noise_norm; nbins = 20)
Теперь попробуем объединить шум и пилообразный сигнал, но сделаем это таким образом, чтобы уровень шума возрастал с уровнем сигнала. Умножим шум на пилообразный сигнал и просуммируем нарастающую "пилу" с нарастающим шумом:
amp_noise = noise .* sawtooth;
noisy_sawtooth = sawtooth .+ (0.25 .* amp_noise);
plot(t, noisy_sawtooth)
Апериодический сигнал и пачка импульсов
В качестве примера непериодического сигнала рассмотрим распространённый sinc-импульс. Мы получаем его при помощи тригонометрической функции и входного вектора значений угла в радианах. В примере мы генерируем импульс в пределах от -2pi до +2pi (с шагом pi/16).
drad = pi/16;
rads = -2*pi:drad:2*pi-drad;
one_sinc = sinc.(rads);
plot(rads, one_sinc, title = "sinc-импульс")
Теперь создадим пачку sinc-импульсов, следующих с определённым периодом. Для этого добавим к вектору одного импульса некоторое количество нулей в конце, и размножим результат функцией repeat. Также добавим нашей пачке импульсов затухание - в этом нам поможет функция LinRange:
zero_padding = [one_sinc; zeros(2*length(one_sinc))]; # sinc-импульс с нулями
four_sinc = repeat(zero_padding, 4); # пачка импульсов без затухания
decay = LinRange(1, 0.4, length(four_sinc)); # сигнал затухания
pulse_train = four_sinc .* decay; # пачка импульсов с затуханием
plot(pulse_train, title = "Пачка sinc-импульсов с затуханием")
Частотно-модулированный сигнал
В заключении рассмотрим частотно-модулированный сигнал с линейной частотной модуляцией (ЛЧМ). Для создания такого сигнала воспользуемся функцией chirp:
LFM = chirp(stoptime, fs, 0, 200; method = "linear");
plot(t, LFM)
Алиасинг при несоблюдении условия теоремы Котельникова
Визуализируем эффект "наложения" отсчётов сигнала с большей основной частотой на отсчёты сигнала с меньшей частотой. Дискретные отсчёты "принадлежат" обоим сигналам, но после дискретизации мы увидим лишь сигнал меньшей частоты:
tvec = 0:0.01:10;
x1 = cos.(pi/5*tvec);
x2 = cos.(11*pi/5*tvec)
plot(tvec,x2, linewidth=3, label = "Исходный сигнал")
tnew = 0:10;
y1 = cos.(11*pi/5*tnew)
plot!(tnew,y1, line=:stem, marker=:circle, linewidth=3, label = "Дискретизация")
plot!(tvec,x1, linewidth=4, label = "Наложение", title = "Наложение, сигнал с частотой в 11 раз меньше")
Если алиасинг уже произошёл, бороться с ним какими-либо методами пост-обработки бесполезно. Необходимо не допускать его возникновения, правильно выбирая частоту дискретизации и используя антиалиасинговые фильтры.