День 3 Летней школы Julia
Часть 4: цифровые фильтры
Цифровые фильтры позволяют подавлять или усиливать отдельные составляющие, или полосы, спектра сигнала.
Рассмотрим пример синтеза КИХ и БИХ-фильтров, а также применение системного объекта фильтра.
# Это нужно выполнить, если не установлены необходимые библиотеки
#Pkg.add("DSP")
#Pkg.add("WAV")
#Pkg.add("Base64")
#Pkg.add("FFTW")
#Pkg.add("FindPeaks1D")
#Pkg.add("Statistics")
using DSP, WAV, Base64, FFTW, FindPeaks1D, Statistics
Исходные данные для обработки
В качестве тестового сигнала мы вновь берём аудиозапись гитарной струны. Но в этот раз мы уже знаем, что это нота Ля малой октавы с основной частотой 220 Гц.
audio_data, fs = wavread("guitar.wav");
nsamples = length(audio_data);
dt = 1/fs;
t = 0:dt:((nsamples - 1) .* dt);
sig = vec(audio_data);
Для прослушивания аудио воспользуемся дополнительной функцией 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
Прослушаем звук колеблющейся струны:
audioplayer(sig,fs)
И отобразим спектр сигнала:
S = DSP.welch_pgram(sig, onesided=true, fs=fs);
plot(freq(S), pow2db.(power(S)))
Выделение первой гармоники
Синтезируем простой эллиптический фильтр нижних частот (ФНЧ), подавляющий частоты выше 360 Гц. Это позволит нам отсечь высшие гармоники. Отфильтруем исходный сигнал и прослушаем результат:
myfilt = digitalfilter(Lowpass(2*360/fs), Elliptic(10, 0.1, 60))
Полезно оценить результат синтеза коэффициентов, визуализировав амплитудно-частотную (АЧХ) и фазо-частотную (ФЧХ) характеристики фильтра, а также его импульсную характеристику:
H, w = freqresp(myfilt);
freq_vec = fs*w/(2*pi);
plot(freq_vec, pow2db.(abs.(H))*2, linewidth=3, title = "АЧХ ФНЧ")
phi, w = phaseresp(myfilt);
plot(freq_vec, rad2deg.(phi/2), linewidth=3, title = "ФЧХ ФНЧ")
myimpresp = impresp(myfilt);
plot(myimpresp, line=:stem, marker=:circle)
Отфильтруем сигнал функцией filt, послушаем результат:
A220 = filt(myfilt, sig);
audioplayer(A220 * 2, fs)
Сравним спектры исходного и отфильтрованного сигнала:
S220 = DSP.welch_pgram(A220, onesided = true, fs = fs);
plot(freq(S), pow2db.(power(S)))
plot!(freq(S220), pow2db.(power(S220)), linewidth = 3)
Отобразим результат фильтрации во временной области:
plot(t, A220)
Теперь попробуем синтезировать полосно-пропускающий КИХ-фильтр, который позволит выделить вторую гармонику в 440 Гц:
ntaps = 100;
b = digitalfilter(Bandpass(2*360/fs, 2*500/fs), FIRWindow(hanning(ntaps)));
myFIR = PolynomialRatio(b,[1]);
H2, w2 = freqresp(myFIR);
plot(freq_vec, pow2db.(abs.(H2))*2, linewidth=3, title = "АЧХ ППФ", ylim=(-140,10))
Результат фильтрации:
A440 = filt(myFIR, sig);
audioplayer(A440, fs)
А также сравнение спектров:
S440 = DSP.welch_pgram(A440, onesided = true, fs = fs);
plot(freq(S), pow2db.(power(S)))
plot!(freq(S440), pow2db.(power(S440)), linewidth=3)
Применение системного объекта фильтра
chunk_len = 800;
nchunks = 10;
sig_crop = sig[1:chunk_len*nchunks]
input_matrix = reshape(sig_crop,(chunk_len,nchunks))
output_matrix = zeros(size(input_matrix));
for i = 1:nchunks
output_matrix[:,i] = filt(myFIR, input_matrix[:,i])
end
out_vector = reshape(output_matrix,chunk_len*nchunks);
t_crop = t[1:length(sig_crop)];
plot(t_crop, out_vector)
audioplayer(out_vector * 4, fs)
Инициализация параметров и применение системного объекта КИХ-фильтра:
fir_SO = EngeeDSP.DescretFIRFilter()
fir_SO.Coefficients = b
output_matrix = zeros(size(input_matrix));
EngeeDSP.setup!(fir_SO, input_matrix[1,:])
for i = 1:nchunks
output_matrix[:,i] = EngeeDSP.step!(fir_SO, input_matrix[:,i])
end
out_vector_SO = reshape(output_matrix,chunk_len*nchunks);
plot(t_crop, out_vector_SO)
audioplayer(out_vector_SO * 4, fs)