День 3 Летней школы Julia
Часть 5: многоскоростная обработка
Многоскоростными называются дискретные системы, содержащие сигналы с разной частотой дискретизации (ЧД). В большинстве случаев рассматриваются системы, изменяющие частоту дискретизации входного сигнала. Её можно повышать и понижать в целое число раз, а также изменять в дробное число раз (m/n, где m и n - целые).
Рассмотрим важность применения фильтров при изменении частоты дискретизации сигнала на примере повышения частоты в 2.5 раза.
# Это нужно выполнить, если не установлены необходимые библиотеки
#Pkg.add("DSP")
#Pkg.add("WAV")
#Pkg.add("Base64")
using DSP, Base64, WAV, DelimitedFiles
Исходные данные для обработки
В качестве тестового сигнала сгенерируем одну синусоиду основной частоты 100 Гц и с частотой дискретизации в 2 кГц. Задача - поднять частоту дискретизации синусоиды в 2.5 раза то есть фактически выполнить операцию интерполяции.
fs = 2000;
f0 = 100;
dt = 1/fs;
stoptime = 0.1;
t = 0:dt:stoptime-dt;
sig = sin.(2*pi*t*f0);
plot(t, sig, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
Убедимся, что на спектре есть только одна значимая компонента на частоте 100 Гц:
S = DSP.periodogram(sig, onesided=true, fs=fs);
plot(freq(S), pow2db.(power(S)))
Коэффициент интерполяции - L = 5
Для начала мы выполняем операцию повышения частоты дискретизации за счёт внесения нулей между отсчётами сигнала:
L = 5;
up = zeros.(length(sig)*L);
up[1:L:end] = sig;
dtL = 1/(fs*L);
tL = 0:dtL:stoptime-dtL;
plot(t, sig, line=:steppost)
plot!(tL, up, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
Форма исходной синусоиды сильно искажена, и также мы можем наблюдать спектральные копии на расширенной первой зоне Найквиста в спектре сигнала:
PL = DSP.periodogram(up, onesided=true, fs=fs*L);
plot(freq(PL), pow2db.(power(PL)), legend = false)
Проредив подобный сигнал, выкинув каждый второй элемент, мы получим требуемую частоту дискретизации, но форма и спектр по-прежнему будут искажены:
M = 2;
down = up[1:M:end];
dtM = 1/(fs*L/M);
tM = 0:dtM:stoptime-dtM;
plot(tM, down, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
Спектральные копии, "завернувшиеся" в первую зону Найквиста:
SM = DSP.periodogram(down, onesided=true, fs=fs*L/M);
plot(freq(SM), pow2db.(power(SM)), legend = false)
Фильтр нижних частот, исполняющий функцию анти-имаджингового при операции повышения частоты дискретизации, а также анти-алиасингового для операции понижения частоты дискретизации:
b = digitalfilter(Lowpass(2*1000/(fs*L)), FIRWindow(hanning(64)));
myfilt = PolynomialRatio(b,[1]);
H, w = freqresp(myfilt);
freq_vec = fs*L*w/(2*pi);
plot(freq_vec, pow2db.(abs.(H))*2, linewidth=3, ylim=(-120,10))
interp = filt(myfilt, up) * L;
plot(tL, interp, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.02))
decim = interp[1:M:end];
plot(t, sig, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
plot!(tM, decim, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
b2 = vec(readdlm("num.txt"));
myfilt2 = PolynomialRatio(b2,[1]);
interp2 = filt(myfilt2, up) * L;
decim2 = interp2[1:M:end];
plot(t, sig, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
plot!(tM, decim2, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
Применение функции resample из состава DSP.jl:
result = resample(sig, L/M);
plot(t, sig, line=:steppost, markersize = 2, xlim=(0,0.01))
plot!(tM[1:end-1], result, line=:steppost, marker=:dot, markersize = 2, xlim=(0,0.01))
Эффект алиасинга при понижении частоты дискретизации:
Для прослушивания аудио воспользуемся дополнительной функцией audioplayer:
function audioplayer(s, fs);
buf = IOBuffer();
wavwrite(s, buf; Fs=fs);
data = base64encode(unsafe_string(pointer(buf.data), buf.size));
markup = """<audio controls="controls" {autoplay}>
<source src="data:audio/wav;base64,$data" type="audio/wav" />
Your browser does not support the audio element.
</audio>"""
display("text/html", markup);
end
Меняя значение частоты дискретизации при помощи слайдера, можно наблюдать изменение основной частоты тона при несоблюдении условия теоремы Котельникова:
fs = 6000 # @param {type:"slider",min:3000,max:6000,step:500}
dt = 1/fs;
t = 0:dt:1;
f0 = 2400;
x = sin.(2*pi*f0*t);
S = DSP.welch_pgram(x, onesided=true, fs=fs);
audioplayer(x,fs)
plot(freq(S), pow2db.(power(S)), legend = false, size=(800,200))