Модуль 4 - Школа анализа инженерных данных

# Вычисление признаков
function вычислить_признаки(сигнал, fs)
    N = length(сигнал)
    сигнал_центр = сигнал .- mean(сигнал)
    abs_сигнал = abs.(сигнал_центр)
    
    # Временные признаки
    rms_val = sqrt(mean(сигнал_центр .^ 2))
    abs_mean = mean(abs_сигнал)
    
    признаки_время = [
        mean(сигнал), # среднее
        median(сигнал), # медиана
        rms_val,        # RMS
        std(сигнал_центр), # СКО
        var(сигнал_центр), # дисперсия
        maximum(сигнал),   # максимум
        minimum(сигнал),   # минимум
        maximum(сигнал) - minimum(сигнал), # размах
        skewness(сигнал_центр), # асимметрия
        kurtosis(сигнал_центр), # эксцесс
        maximum(сигнал) / (rms_val + 1e-10), # пик-фактор
        maximum(сигнал) / (abs_mean + 1e-10), # импульс-фактор
        rms_val / (abs_mean + 1e-10), # форма-фактор
        sum(сигнал_центр .^ 2), # энергия
    ]
    
    # Спектральные признаки
    спектр = abs.(fft(сигнал_центр))
    одност = спектр[1:N÷2]
    частоты = (0:N÷2-1) * fs / N
    норм = одност / (sum(одност) + 1e-10)
    кум_спектр = cumsum(норм)
    макс_спектра = maximum(одност)
    домин_частота = частоты[argmax(одност)]
    спектр_центр = sum(частоты .* норм)
    спектр_ширина = sqrt(sum(((частоты .- спектр_центр).^2) .* норм))
    энтропия = -sum(норм[норм .> 0] .* log.(норм[норм .> 0]))
    геом_среднее = exp(mean(log.(норм[норм .> 0])))
    спектр_плоскостность = геом_среднее / (mean(норм) + 1e-10)
    средняя_частота = sum(частоты .* норм)
    idx_медиана = findfirst(>=(0.5), кум_спектр)
    медианная_частота = idx_медиана === nothing ? 0.0 : частоты[idx_медиана]
    idx_rolloff = findfirst(>=(0.85), кум_спектр)
    роллоф = idx_rolloff === nothing ? 0.0 : частоты[idx_rolloff]
    спектр_энергия = sum(одност .^ 2)
    спектр_ско = std(одност)
    спектр_асимм = skewness(одност)
    спектр_эксцесс = kurtosis(одност)
    
    признаки_спектр = [
        макс_спектра, # максимум спектра
        домин_частота, # доминирующая частота
        спектр_центр,  # спектральный центр
        спектр_ширина, # спектральная ширина
        энтропия,      # спектральная энтропия
        спектр_плоскостность, # спектральная плоскостность
        средняя_частота,      # средняя частота
        медианная_частота,  # медианная частота
        роллоф,          # spectral roll-off 85%
        спектр_энергия,  # спектральная энергия
        спектр_ско,      # спектральное СКО
        спектр_асимм,    # спектральная асимметрия
        спектр_эксцесс,  # спектральный эксцесс
    ]
    
    # Энергия частотных полос
    полосы = [
        (0,500), (500,1000), (1000,1500), (1500,2000), (2000,2500),
        (2500,3000), (3000,3500), (3500,4000), (4000,4500), (4500,5000)
    ]
    энергии_полос = [sum(одност[(частоты .>= н) .& (частоты .< в)] .^ 2) for (н, в) in полосы]
    
    # Признаки огибающей
    огиб = abs.(DSP.hilbert(сигнал_центр))
    спектр_огиб = abs.(fft(огиб .- mean(огиб)))[1:N÷2]
    норм_огиб = спектр_огиб / (sum(спектр_огиб) + 1e-10)
    
    макс_огиб = maximum(спектр_огиб)
    домин_огиб = частоты[argmax(спектр_огиб)]
    энергия_огиб = sum(спектр_огиб .^ 2)
    энтропия_огиб = -sum(норм_огиб[норм_огиб .> 0] .* log.(норм_огиб[норм_огиб .> 0]))
    центр_огиб = sum(частоты .* норм_огиб)
    
    признаки_огиб = [
        макс_огиб,     # максимум спектра огибающей
        домин_огиб,    # доминирующая частота огибающей
        энергия_огиб,  # энергия спектра огибающей
        энтропия_огиб, # энтропия огибающей
        центр_огиб,    # спектральный центр огибающей
    ]
    
    # Дополнительные признаки для ролика
    # Статистики огибающей во времени
    огиб_среднее = mean(огиб)
    огиб_ско = std(огиб)
    огиб_макс = maximum(огиб)
    огиб_эксцесс = kurtosis(огиб)
    
    # Пересечения нуля
    пересечения_нуля = sum(diff(sign.(сигнал_центр)) .!= 0) / N
    
    # Энергетические отношения
    вч_энергия = sum(энергии_полос[8:10])
    сч_энергия = sum(энергии_полос[4:7])
    нч_энергия = sum(энергии_полос[1:3]) + 1e-10
    
    отн_вч_нч = вч_энергия / нч_энергия
    отн_огиб_спектр = энергия_огиб / (спектр_энергия + 1e-10)
    импульсность = maximum(abs_сигнал) / (rms_val + 1e-10)
    
    # Анализ амплитудного распределения
    сортировка = sort(abs_сигнал)
    p25 = сортировка[div(N, 4)]
    p50 = сортировка[div(N, 2)]
    p75 = сортировка[div(3*N, 4)]
    
    малые = abs_сигнал[abs_сигнал .< p25]
    большие = abs_сигнал[abs_сигнал .> p75]
    
    доля_малых = length(малые) / N
    отн_малых_больших = mean(малые) / (mean(большие) + 1e-10)
    отн_энергии_малых = sum(малые .^ 2) / (sum(большие .^ 2) + 1e-10)
    отн_p75_p25 = p75 / (p25 + 1e-10)
    
    return vcat(
        признаки_время,    
        признаки_спектр,       
        энергии_полос,          
        признаки_огиб,         
        [огиб_среднее, огиб_ско, огиб_макс, огиб_эксцесс,
         пересечения_нуля, отн_вч_нч, отн_огиб_спектр, импульсность,
         доля_малых, отн_малых_больших, отн_энергии_малых, отн_p75_p25]
    )
end