• Регистрация
masteroforest
masteroforest 0.00
н/д

Преобразование Фурье в MATLAB

21.12.2020

Подскажите, пожалуйста, относительно Фурье преобразования в среде MATLAB. Сам несколько дней пытаюсь понять через гугл и периодические журналы... Есть набор ЭЭГ-сигналов: 500 значений за 2 сек (частота 250 Гц). Мне надо просто разложить их в спектр (насколько понимаю от 0 до 125) и получить значения 2-х полос: от 0.5 до 14 Гц и от 14 до 40 Гц. Но даже разложение в спектр, полагаю с ошибкой:

% init
fs = 250;

% load
da = da1_all;
x = da(1:256,1);

% Построение графика спектральной плотности
fx = fft(x);
sfx = fftshift(fx);
f = linspace(-fs/2,fs/2,length(x));
area(f,abs(sfx));
title('Амплитудный спектр сигнала');
xlabel('Частота(Гц)');
ylabel('Амплитуда');

Теги

    21.12.2020

    Ответы

    • Andrey Ermakov
      Andrey Ermakov +56.93
      21.12.2020 12:14

      Может быть это поможет

      https://www.mathworks.com/matlabcentral/fileexchange/50082-eeg-data-plotting-power-spectrum-spectrogram-frequency-spectrum-of-alpha-beta-delta-and-theta?s_tid=srchtitle

      • masteroforest
        masteroforest0.00
        21.12.2020 13:31

        Благодарю. Но там самого кода преобразования в спектр не нашёл(фигурирует только fs/2).

      • aBoomest
        aBoomest+942.89
        21.12.2020 12:19

        1. sfx = fftshift(fx); - это зачем?
        2. Расскажите хоть что-то о своем сигнале (форма, гармоники, и т.д.), чтоб можно было анализировать? Или там какие-то сейсмические данные о которых вы ничего не знаете? А то, то что у вас на картинке, от этого ни жарко, ни холодно, ибо информации входной нет. Ну или код выложите, при условии если он не огромных размеров :).

        • masteroforest
          masteroforest0.00
          21.12.2020 13:46

          Для примера (тут 8 каналов). Спектральное представление нужно для каждого индивидуально. Также здесь 2 сек, а сейчас работаю с 4 сек. временными рядами. 250 отсчётов в 1 сек. По-хорошему, целевой сигнал находится в диапазоне частот 0,5 - 40 Гц (вне диапазона шумы электросети и др...). fftshit, как понял из периодики, нужен для сдвига нулевой частоты в начало координат (т. е. результат не был симметричным).

          P.S. Да, сигнал я не синтезирую, а снимаю со своей головы (т.е. кода никакого нет).

          • aBoomest
            aBoomest+942.89
            21.12.2020 15:21

            fftshift думаю не нужен.
            сохраните 1 сигнал в mat файле и выложите , проверить. или эксель, или в чем это у вас.

            ps: а вобще в хелпе на fft всю жизнь был толковый пример.

            psps: с головы? )))))))))

        • masteroforest
          masteroforest0.00
          21.12.2020 19:06

          С поверхности головы, точнее не скажешь)

          Приложил сигнал по одному каналу (1 из 8 как с графика выше).

          А где в хелпе? Пробовал по примеру с офсайта https://se.mathworks.com/help/matlab/ref/fft.html

          • aBoomest
            aBoomest+942.89
            22.12.2020 09:13

            Так про него и говорю.

             

            • aBoomest
              aBoomest+942.89
              22.12.2020 09:53

              Вы наблюдаете ВЧ составляющие в вашем сигнале? Вот они и вылезают. 

              Не ясно следующее:
              1. Вы прислали самый не показательный сигнал (нижнюю кривую) - она самая не изменяющаяся.
              2. Увас на рисунке по оси ОХ написано - отсчеты. При этом максимальное значение 2е4. Т.е. 20000 - отсчетов. А в файле, что выложили, только 500. 

          • masteroforest
            masteroforest0.00
            22.12.2020 13:05

            Сигналы с графика и с mat-файла, действительно разные. Они сделаныв разное время, но при одинаковых условиях. Для графика нужна была показательность, чего нельзя сделать на коротком интервале. Задача перевода в спектр встала несколько позднее, тогда начал деалть отсчётыпо 2 сек. А прислал я не нижнюю кривую, а кривую 1 (наоборот, то ли верхнюю, то ли из центра).

            Высокочастотные, да. Имеющий смысл диапазон находится до 40 Гц. Выше 50 Гц начинаются наводки электросети (поэтому собирался обрезать).

            Вообще, спасибо за наглядные графики. Можете сориентировать, что скорректировать в моем коде(в идеале, конечно, ваш пример кода, по которому нижний график получился)?

            • aBoomest
              aBoomest+942.89
              22.12.2020 16:33

              А прислал я не нижнюю кривую, а кривую 1 (наоборот, то ли верхнюю, то ли из центра).

              Но у вас там есть "волны", а на рис.предыдущего поста - просто "прямая с шумом".

              UPD: Быть не может, у вас там только одна кривая с отрицательными значениями.

              Высокочастотные, да. Имеющий смысл диапазон находится до 40 Гц. Выше 50 Гц начинаются наводки электросети (поэтому собирался обрезать).
              - сигнал
              - полный спектр
              - спектр 0-40 Гц

              C1 - из вашего файла.
              Вот мне не жалко, честно, но кроме
              - чтения вашего файла,
              - функции Y = fft(C1) ,
              - затем вывода графика,

              . . . ничего больше нет. Никаких других хитрых функций и замысловатостей, которые у вас присутствуют. Попробуйте все лишнее выкинуть. Заодно Вам же понятнее станет. Если уж никак не получится, то конечно выложу. :)
              PS: Но надо конечно "фильтрануть" всё лишнее, что вам не интересно, постояннное смещение особенно, иначе на фоне вы ничего полезного не увидите. Хотя я не знаю, если вам постоянка нужна, то графики, что нарисовал, не верные.

              • masteroforest
                masteroforest0.00
                24.12.2020 15:02

                Добрый день!

                Да, пример кода всё-таки не помешает. Попробовал упростить, вот что получилось:

                Про постоянную составляющую: предполагаю, она всё-таки не нужна. С фильтрацией и удалением шумов, конечно, ещё предстоит работа...

                Не знаю точно, почему на некоторых каналах разность потенциалов (амплитуда) отрицательная. Не стал углубляться, т.к. счёл ключевыми параметрыми свойства волны. Ещё раз оговорюсь, сигнал, котрый отображен на верхней картинке снимался ранее (был программынй сбой и проходили только 8 каналов). Насколько помню, его протяженность 80 сек. Сигнал из mat-файла по 1 из 16 каналов, с которыми работаю сейчас. Но на коротких промежутках он действительно выглядит как прямая с шумом:

                • aBoomest
                  aBoomest+942.89
                  24.12.2020 20:46

                  Код у вас наверняка правильный, судя по тычку на 200 Гц, там ошибиться сложно. Удалите сперва постоянное смещение из сигналов (безо всяких фильтров). Самый простой деревенский способ:

                  Signal = Signal / mean(Signal);

                  Думаю сразу что-то да прояснится в спектре.