266 lines
8.8 KiB
Julia
266 lines
8.8 KiB
Julia
include("./zaf.jl")
|
||
using .zaf
|
||
using WAV
|
||
using Statistics
|
||
using Plots
|
||
using Tables
|
||
using CSV
|
||
using DSP
|
||
using SignalAnalysis
|
||
using DataFrames
|
||
using LinearAlgebra
|
||
using Logging
|
||
|
||
global file_frame = ""
|
||
|
||
global savedata = false #Сохранять ли при обработке спектограммы
|
||
|
||
detect = true #Флаг необходимости сравнения с эталонами
|
||
|
||
# Настройки CQT преобразования
|
||
|
||
octave_resolution = 24
|
||
minimum_frequency = 55
|
||
maximum_frequency = 9520
|
||
time_resolution = 25
|
||
|
||
# Размеры кадра в сек
|
||
|
||
window_duration = 0.04;
|
||
|
||
# Уровень энергии для порогового фильтра
|
||
|
||
energy_level = 50
|
||
|
||
# Количество кадров для сглаживания скользящей средней слепка аудио-данных
|
||
|
||
m_offset = 20
|
||
|
||
function err(mes)
|
||
redirect_stdio(stderr="audio.log") do
|
||
logger = ConsoleLogger(stderr)
|
||
@error mes
|
||
with_logger(logger) do
|
||
println(mes)
|
||
end
|
||
end
|
||
exit(2)
|
||
end
|
||
|
||
# Загрузка сохраненного файла данных очищенных окон с преобразованием в среднеквадратичное отклонение по частотным полосам
|
||
|
||
function load(filename)::Vector{Any}
|
||
if !isempty(filename)
|
||
try
|
||
data = CSV.read(filename, DataFrame, skipto=2)
|
||
return data[!, :Column1]
|
||
catch e
|
||
println("Ошибка чтения csv-файла")
|
||
end
|
||
else
|
||
err("Не указан входной файл")
|
||
end
|
||
end
|
||
|
||
function sets()::Dict{String, Vector{Float64}}
|
||
datasets = Dict{String, Vector{Float64}}()
|
||
for file in readdir("datasets")
|
||
datasets[replace(file,".csv" => "")] = load("datasets/"*file)
|
||
end
|
||
return datasets
|
||
end
|
||
|
||
# Сравнение частотного снимка с сохраненными эталонами
|
||
|
||
function determ(data)
|
||
for (name,frame) in sets()
|
||
if cor(data[70:end], frame[70:end]) < 0.99 #Если корреляция по Пирсону по ВЧ низкая, то сигнал невалидный
|
||
continue
|
||
end
|
||
dataLF = normalize(data[1:70])
|
||
frame = normalize(frame[1:70])
|
||
if sum(broadcast(abs, dataLF-frame)) < 0.2 #Если сумма асболютных отклнений по НЧ низкая, то фрейм совпадает
|
||
return name
|
||
end
|
||
end
|
||
return
|
||
end
|
||
|
||
# Скользящее среднее
|
||
|
||
function m_average(arr::Vector{Float64}, n::Number)::Vector{Any}
|
||
so_far = sum(arr[1:n])
|
||
out = zero(arr[n:end])
|
||
out[1] = so_far
|
||
for (i, (start, stop)) in enumerate(zip(arr, arr[n+1:end]))
|
||
so_far += stop - start
|
||
out[i+1] = so_far
|
||
end
|
||
return out
|
||
end
|
||
|
||
# Нормализация уровня
|
||
|
||
function normalizeWAV(audio_input::Matrix{Float64},freq::Number)
|
||
return mean(audio_input, dims=2), freq
|
||
end
|
||
|
||
# Разбитие на окна
|
||
|
||
function windows(freq::Number)
|
||
local window_length = nextpow(2, ceil(Int, window_duration*freq))
|
||
local window_function = zaf.hamming(window_length, "periodic") # оконная аналитическая функция
|
||
return window_length, convert(Int, window_length/2), window_function
|
||
end
|
||
|
||
# Фильтрация окон по уровню минимальной энергии звука в каждом окне
|
||
|
||
function filter(audio_input::Matrix{Float64},window_length::Number)::Matrix{Float64}
|
||
output_signal = zeros(1,1)
|
||
number_samples = length(audio_input)
|
||
padding_length = floor(Int, window_length / 2)
|
||
number_times =
|
||
ceil(
|
||
Int,
|
||
((number_samples + 2 * padding_length) - window_length) /
|
||
step_length,
|
||
) + 1
|
||
local i = 0
|
||
for j = 1:(number_times-10) # 10 последних кадров отбарсываем
|
||
frame = audio_input[i+1:i+window_length]
|
||
if energy(frame) > energy_level
|
||
output_signal = [output_signal;frame]
|
||
end
|
||
i = i + step_length
|
||
|
||
end
|
||
return output_signal
|
||
end
|
||
|
||
# Преобразование Фурье
|
||
|
||
function stft(audio::Matrix{Float64},freq::Number)
|
||
local window_length,step_length,window_function = windows(freq)
|
||
return zaf.stft(audio, window_function, step_length), window_length, step_length, window_function
|
||
end
|
||
|
||
# Спектограмма разложения ряда Фурье
|
||
|
||
function spectrogram(audio::Matrix{Float64},freq::Number)::Matrix{Float64}
|
||
local audio_stft, window_length = stft(audio,freq)
|
||
local audio_spectrogram = abs.(audio_stft[2:convert(Int, window_length/2)+1, :]) # Отбрасываем постоянную составляющую
|
||
local number_samples = length(audio_input)
|
||
local xtick_step = 1
|
||
local ytick_step = 1000
|
||
local plot_object = zaf.specshow(audio_spectrogram, number_samples, freq, xtick_step, ytick_step)
|
||
heatmap!(title = "Spectrogram (dB)", size = (990, 600))
|
||
savefig(plot_object, file_frame*".spectr.png")
|
||
CSV.write(file_frame*".spectr.csv", Tables.table(audio_spectrogram))
|
||
return audio_spectrogram
|
||
end
|
||
|
||
# CQT-октавное вевлет-преобразование
|
||
|
||
function cqt(audio::Matrix{Float64},freq::Number)::Matrix{Float64}
|
||
local cqt_kernel = zaf.cqtkernel(freq, octave_resolution, minimum_frequency, maximum_frequency)
|
||
cqt_spectrogram = zaf.cqtspectrogram(audio, freq, time_resolution, cqt_kernel)
|
||
if savedata
|
||
CSV.write(file_frame*".cqt.csv", Tables.table(cqt_spectrogram))
|
||
local xtick_step = 1
|
||
local plot_object = zaf.cqtspecshow(cqt_spectrogram, time_resolution, octave_resolution, minimum_frequency, xtick_step)
|
||
heatmap!(title = "CQT spectrogram (dB)", size = (990, 600))
|
||
savefig(plot_object, file_frame*".cqt.png")
|
||
end
|
||
return cqt_spectrogram
|
||
end
|
||
|
||
# Преобразование вевлет-спектрограммы в спектр среднего значения по октавам со сглаживанием
|
||
|
||
function raft(spectr::Matrix{Float64})::Vector{Any}
|
||
freqs, frames = size(spectr)
|
||
raft_data = zeros(freqs)
|
||
for i in 1:freqs
|
||
raft_data[i] = mean(spectr[i, :])
|
||
end
|
||
return m_average(raft_data,m_offset)
|
||
end
|
||
|
||
# Cохранение образа записи
|
||
|
||
function save(audio_data::Vector{Any},freq::Number)::Vector{Any}
|
||
CSV.write(file_frame*".csv", Tables.table(audio_data))
|
||
numbers = length(audio_data)
|
||
xtick_locations = [0:octave_resolution:numbers;]
|
||
xtick_labels = convert(
|
||
Array{Int},
|
||
minimum_frequency * 2 .^ (xtick_locations / octave_resolution),
|
||
)
|
||
plot_object = plot(audio_data, xticks = (xtick_locations, xtick_labels))
|
||
plot!(title = "Sound form", size = (990, 600))
|
||
savefig(plot_object, file_frame*".png")
|
||
return audio_data
|
||
end
|
||
|
||
|
||
# Считываем данные из WAV файла и нормализуем их, частота дискретизации получается из заголовка WAV-файла
|
||
|
||
if !isempty(ARGS)
|
||
global file_frame = replace(ARGS[1],".wav" => "")
|
||
if length(ARGS)>1 && ARGS[2] == "-s"
|
||
savedata = true
|
||
end
|
||
if length(ARGS)>1 && ARGS[2] == "-d"
|
||
println("Режим формирования нового эталонного датасета")
|
||
if ARGS[3] isa String
|
||
savedata = true
|
||
global file_frame = "datasets/"*ARGS[3]
|
||
detect = false
|
||
else
|
||
err("Не указано название датасета")
|
||
end
|
||
end
|
||
else
|
||
err("Не указан входной файл")
|
||
end
|
||
|
||
audio_input, freq = wavread(ARGS[1])
|
||
|
||
audio_input, freq = normalizeWAV(audio_input, freq)
|
||
|
||
# Исходя из частоты дискретизации получаем длину кадров и шаг между ними
|
||
|
||
window_length,step_length,window_function = windows(freq)
|
||
|
||
# Фильтруем кадры, отбрасывая кадры с низкой энергией
|
||
|
||
audio = filter(audio_input,window_length)
|
||
|
||
frames, freqs = size(audio)
|
||
|
||
if frames<5
|
||
err("Недостаточная энергия сигнала: "*ARGS[1])
|
||
end
|
||
|
||
# Осущетсвляем CQT-преобразование и сохранение спектограммы и массива данных
|
||
|
||
cqt_spectrogram = cqt(audio,freq)
|
||
|
||
# Преобразуем спектограмму в с спектр среднеквадратичных отклонений
|
||
|
||
raft_data = raft(cqt_spectrogram)
|
||
|
||
# Cохраняем результат
|
||
|
||
if savedata
|
||
#spectrogram(audio_input,freq)
|
||
save(raft_data,freq)
|
||
end
|
||
|
||
# Сравниваем с эталонами и выводим результат в stdout
|
||
|
||
if detect
|
||
println(determ(raft_data))
|
||
end
|
||
|
||
|