Выявление отклонений от нормы параметров мозговых волн человека

                           2,     0,     2 ...                         

 ...%                     10     11     12                          

                           0,     2,     2 ...                                                    

 ...%                     13     14     15                          

                           2,     2,     2 ...                         

 ...%                     16     17     18                          

                           2,     0,     0         ];

                    

false_negatives = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];

 

fs = 100; % частота дискретизации

 

for i=1:max(size(healthy))

    'Healthy patient: '

    i

    data = dlmread(healthy{i, 1});

    for j=1:18

        if (channel_indicators(1,j) > 0)

            period = max(size(data));

            n = round(20/fs*period);

 

            ffc = filtered_fourier_coeffs( data(:,j), n, fs, 1, period, 0.5 );

            ffc = [ffc(:,1)/2/pi, sqrt(ffc(:,2).*ffc(:,2)+ffc(:,3).*ffc(:,3))];

 

            lower = lower_boundaries(1, j);

            left = left_boundaries(1, j);

            right = right_boundaries(1, j);

            max_num_of_hits = max_nums_of_hits(1, j);

           

            ffc = ffc(round(left*period/fs+1):round(right*period/fs-1), :);           

            num_of_hits = length(find(ffc(:,2) > lower));

            if (num_of_hits > max_num_of_hits)

                false_negatives(1,j) = false_negatives(1,j) + 1;

            end;

           

        end;

    end;       

end;

 

false_positives = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];

for i=1:max(size(sick))

    'Sick patient: '

    i

    data = dlmread(sick{i, 1});

    for j=1:18

        if (channel_indicators(1,j) > 0)

            period = max(size(data));

            n = round(20/fs*period);

 

            ffc = filtered_fourier_coeffs( data(:,j), n, fs, 1, period, 0.5 );

            ffc = [ffc(:,1)/2/pi, sqrt(ffc(:,2).*ffc(:,2)+ffc(:,3).*ffc(:,3))];

 

            lower = lower_boundaries(1, j);

            left = left_boundaries(1, j);

            right = right_boundaries(1, j);

            max_num_of_hits = max_nums_of_hits(1, j);

           

            ffc = ffc(round(left*period/fs+1):round(right*period/fs-1), :);

            hold on;

            scatter(ffc(:,1), ffc(:,2), 15, 'red', 'filled');

            num_of_hits = length(find(ffc(:,2) > lower));

            if (num_of_hits <= max_num_of_hits)

                false_positives(1,j) = false_positives(1,j) + 1;

            end;

           

        end;

    end;       

end;

 

false_negatives / max(size(healthy)) * 100

false_positives / max(size(sick)) * 100

 

function y = filtered_fourier_coeffs( func, n, fs, start_time, period, level )

% Возвращает оценку коэффициентов  Фурье для функции,

% заданной рядом значений func в дискретные моменты времени  time.

% Гармоники, для которых  a_n^2+b_n^2 <= level^2, зануляются

%

% fs - частота дискретизации

% period - период сигнала (в  отсчетах)

% func (столбец) - значения функции  в дискретные моменты времени

% n     - номер гармоники,  до которой следует считать  коэффициенты Фурье

% level - пороговое значение  амплитуды гармоник

%

% Возвращаемое значение: [omega*k, a_k, b_k],

% где omega*k - столбец из  частот,

% a_k - столбец из коэффициентов  для cos,

% b_k - столбец из коэффициентов  для sin (k = 0..n)

%

% size(y) - n+1 x 3 (для отфильрованных  гармоник a_n = b_n = 0)

%

% func(t) = a_0/2 + a_1*cos(omega*1*t) + b_1*sin(omega*1*t) + ...

%

    y = fourier_coeffs(func,n, fs, start_time, period);

    filter = y(:,2).*y(:,2)+y(:,3).*y(:,3) <= level*level;

    y([filter ~= filter, filter, filter]) = 0;

end

 

 

function y = fourier_coeffs( func, n, fs, start_time, period )

% Возвращает оценку коэффициентов  Фурье для функции,

% заданной рядом значений func в дискретные моменты времени  time

%

% func (столбец) - значения функции  в дискретные моменты времени  time

% n    - номер гармоники,  до которой следует считать  коэффициенты Фурье

% fs   - частота дискретизации

% period - период (в отсчетах)

%

% Возвращаемое значение: [omega*k, a_k, b_k],

% где omega*k - столбец из  частот,

% a_k - столбец из коэффициентов  для cos,

% b_k - столбец из коэффициентов  для sin (k = 0..n)

%

% size(y) - n+1 x 3

%

% func(t) = a_0/2 + a_1*cos(omega*1*t) + b_1*sin(omega*1*t) + ...

%

l = period/fs/2; % T/2

 

time=(0:1:period-1)'*1/fs;

 

func = func(start_time:start_time+period-1);

 

y = zeros([n+1, 3]);

k = 0;

    while(k <= n)

        y(k+1,1:3)=[k*pi/l, trapz(time, func.*cos(k*pi/l.*time))./l, trapz(time, func.*sin(k*pi/l.*time))./l];

 

        k = k + 1;

    end;

end

 

 

Листинг 2. Программа для построения графиков спектральных амплитуд, полученных из ЭЭГ больных и здоровых людей

healthy={'data/healthy/ASHIN_data.txt',

'data/healthy/Fon-atlas_data.txt',

'data/healthy/Fon-avin_data.txt',

'data/healthy/fon-bajenov_data.txt',

'data/healthy/Fon-balaban_data.txt',

'data/healthy/Gavrilov_data.txt'};

 

sick={'data/sick/parkinson/CHERCHES_data.txt',

'data/sick/parkinson/ESMERZUE_data.txt',

'data/sick/parkinson/IABLONOV_data.txt',

'data/sick/parkinson/KUZNECOV_data.txt',

'data/sick/parkinson/LUCHNIKO_data.txt',

'data/sick/parkinson/PANSHIN_data.txt'};

 

for i=1:18       

    clf;

   

    for j=1:max(size(sick))

        data = dlmread(sick{j,1}, '\t');       

        for k=1:1

            hold on;

            fft_convampl_plot(data((k-1)*4096+1:k*4096,i), 4096, 100, 8, 'red', 1, 20);

        end;

    end;

   

    for j=1:max(size(healthy))

        data = dlmread(healthy{j,1}, '\t');       

        for k=1:1

            fft_convampl_plot(data((k-1)*4096+1:k*4096,i), 4096, 100, 8, 'blue', 1, 20);

        end;

    end;

        

    grid on;

    xlim([1,20]);

    saveas(gcf, [num2str(i) '.jpg']);   

       

end;

hold off;

 

function y = plot_ffc2( data, fs, start_time, period, num_of_harmonics, level, color)

% data - столбец данных

% fs - частота дискретизации

% start_time - отсчет, с которого  будут анализироваться данные

% period - период (в отсчетах) сигнала

% num_of_harmonics - число гармоник, которое следует оставить в  каждом из 

%                    частотных диапазонов

 

    n   = round(20/fs*period+1); % считаем коэфф Фурье для гармоник до 20 Гц

    tStart = tic;

    ['start filtered_fourier_coeffs ...']

    ffc = filtered_fourier_coeffs( data, n, fs, start_time, period, level );

    ['end filtered_fourier_coeefs, ex time = ' toc(tStart)]

   

    ffc = [ffc(:,1)/2/pi, sqrt(ffc(:,2).*ffc(:,2)+ffc(:,3).*ffc(:,3))];

   

    delta = ffc(round(0.5*period/fs):round(4*period/fs), :);

    theta = ffc(round(4*period/fs):round(8*period/fs), :);

    alpha = ffc(round(8*period/fs):round(13*period/fs), :);

    beta = ffc(round(13*period/fs):round(20*period/fs), :);

   

    delta = process_range(delta, num_of_harmonics);

    theta = process_range(theta, num_of_harmonics);

    alpha = process_range(alpha, num_of_harmonics);

    beta  = process_range(beta, num_of_harmonics);

   

    hold on;

    scatter(delta(:,1), delta(:,2), 15, color, 'filled');

    scatter(theta(:,1), theta(:,2), 15, color, 'filled');

    scatter(alpha(:,1), alpha(:,2), 15, color, 'filled');

    scatter(beta(:,1), beta(:,2), 15, color, 'filled');

    hold off;

   

    y = ffc;

end

 

function y = process_range(range, num_of_harmonics)

    [S, I] = sort(range(:,2), 'descend');

    range = [range(I, 1), range(I, 2)];

       

    if (num_of_harmonics > 0 & max(size(range)) > num_of_harmonics)  

        range = [range(1:num_of_harmonics, 1), range(1:num_of_harmonics, 2)];

    end;

   

    for k=1:size(range(:,1))

        if (range(k,2) == 0)

            break;

        end;

    end;

    k = max([0,k-1]);

    y = range(1:k, :);

end

 

Листинг 3. Программа для нахождения матриц ошибок первого и второго рода для метода с использованием функций правдоподобия

healthy={'data/healthy/ASHIN_data.txt'

'data/healthy/Fon-atlas_data.txt'

'data/healthy/Fon-avin_data.txt'

'data/healthy/fon-bajenov_data.txt'

'data/healthy/Fon-balaban_data.txt'

'data/healthy/Gavrilov_data.txt'

'data/healthy/Voronov_data.txt'

};

 

sick={'data/sick/parkinson/PANSHIN_data.txt'

'data/sick/parkinson/ESMERZUE_data.txt'

'data/sick/parkinson/IABLONOV_data.txt'

'data/sick/parkinson/KUZNECOV_data.txt'

'data/sick/parkinson/LUCHNIKO_data.txt'

};

 

sick_t={'data/sick/parkinson/CHERCHES_data.txt'

'data/sick/parkinson/ESMERZUE_data.txt'

'data/sick/parkinson/IABLONOV_data.txt'

'data/sick/parkinson/KUZNECOV_data.txt'

'data/sick/parkinson/LUCHNIKO_data.txt'

'data/sick/parkinson/PANSHIN_data.txt'};

 

healthy_t={'data/healthy/ASHIN_data.txt'

'data/healthy/Fon-atlas_data.txt'

'data/healthy/Fon-avin_data.txt'

'data/healthy/fon-bajenov_data.txt'

'data/healthy/Fon-balaban_data.txt'

'data/healthy/Gavrilov_data.txt'

'data/healthy/Gliken_data.txt'

'data/healthy/Gorchakov-long_data.txt'

'data/healthy/Gorchakov1_data.txt'

'data/healthy/Gorchakov_data.txt'

'data/healthy/Grigorev_data.txt'

'data/healthy/Gunin_data.txt'

'data/healthy/Gusev_data.txt'

'data/healthy/Haritonova_data.txt'

'data/healthy/IVANOV1_data.txt'

'data/healthy/IVANOV2_data.txt'

'data/healthy/Ivanova_data.txt'

'data/healthy/Ivanov_data.txt'

'data/healthy/Jashin_data.txt'

'data/healthy/KILDUSHK_data.txt'

'data/healthy/Kiseleva_data.txt'

'data/healthy/Kozlova_data.txt'

'data/healthy/Kozlov_data.txt'

'data/healthy/Krlehsoev_data.txt'

'data/healthy/Lamkin_data.txt'

'data/healthy/Leonon_data.txt'

'data/healthy/Lilna_data.txt'

'data/healthy/Midzjanovskaja_data.txt'

'data/healthy/Ogurtcova-b_data.txt'

'data/healthy/Ogurtcova_data.txt'

'data/healthy/Ogurtcva-a_data.txt'

'data/healthy/Olferev_data.txt'

'data/healthy/OSTROVS1_data.txt'

'data/healthy/OSTROVS2_data.txt'

'data/healthy/Sazonov1_data.txt'

'data/healthy/Sazonov2_data.txt'

'data/healthy/Sazonov3_data.txt'

'data/healthy/Shonina_data.txt'

'data/healthy/Silaeva1_data.txt'

'data/healthy/Silaeva2_data.txt'

'data/healthy/Stepanuk1_data.txt'

'data/healthy/Stepanuk2_data.txt'

'data/healthy/Stepanuk3_data.txt'

'data/healthy/Voronov_data.txt'

};

 

n = 3000;

false_positives = zeros(max(size(sick)), max(size(healthy)));

false_negatives = zeros(max(size(sick)), max(size(healthy)));

size(false_positives)

 

for i=1:max(size(healthy))

    hdata = dlmread(healthy{i, 1});   

    hdata = norm_data(hdata(1:n, 1:18));

    hcov = edfcov(hdata);

   

    for j=1:max(size(sick))

        sdata = dlmread(sick{j, 1});

        sdata = norm_data(sdata(1:n, 1:18));

        scov = edfcov(sdata);      

       

        'i, j'

        i

        j       

        for k=1:max(size(sick_t))

            tsdata = dlmread(sick_t{k, 1});

            tsdata = norm_data(tsdata(1:n, 1:18));

           

            if (lkfunc(hcov, tsdata) > lkfunc(scov, tsdata))

                false_positives(j,i) = false_positives(j,i)+1;

            end;

        end;

       

        for k=1:max(size(healthy_t))

            thdata = dlmread(healthy_t{k, 1});

            thdata = norm_data(thdata(1:n, 1:18));

           

            if (lkfunc(hcov, thdata) < lkfunc(scov, thdata))

                false_negatives(j,i) = false_negatives(j,i)+1;

            end;

        end;

                       

    end;

end;

 

false_negatives / max(size(healthy_t)) * 100

false_positives / max(size(sick_t)) * 100

 

function y = lkfunc( cov, ndata )

    n = max(size(ndata(:,1)));

    covDet = det(cov);

    iCov = inv(cov);

    S = 0;

    for k=1:n

        S = S - 1/2*ndata(k,:)*iCov*(ndata(k,:)');

    end;

   

    Z = S-9*n*log(2*pi)-n/2*log(covDet);

    %L = exp(S)/sqrt((2*pi)^(18)*covDet)^n;

    y = Z;

end

 

 

function y = edfcov( ndata ) 

    n = max(size(ndata(:, 1)));

    num_of_sensors = max(size(ndata(1, :)));

       

    cov = zeros(num_of_sensors, num_of_sensors);

    for i=1:num_of_sensors

        for j=1:num_of_sensors

            c = 0;

            for k=1:n               

                c = c + ndata(k, i)*ndata(k,j);

            end;

           

            cov(i,j) = c / n;

        end;

    end;

   

    y = cov;

       

end

 

 

function y = norm_data( data )

    n = max(size(data(:,1)));

    num_of_sensors = max(size(data(1, :)));   

   

    for i=1:num_of_sensors

        iMean = 0;

        iVar = 0;

       

        for k=1:n

            iMean = iMean + data(k, i);

            iVar = iVar + data(k, i) * data(k, i);

        end;

       

        iMean = iMean / n;

        iVar = iVar / n - iMean*iMean;

       

        for k=1:n

            data(k, i) = (data(k, i) - iMean) / sqrt(iVar);

        end;

    end;

       

    y = data;

end

 

 

 

1