MATLAB и быстрое преобразование Фурье

По работе неоднократно сталкивался с необходимостью быстро определить наличие в сигнале гармонических составляющих. Часто для примерной оценки достаточно воспользоваться алгоритмом быстрого преобразования Фурье. Тем более, что его реализации есть практически во всех математических пакетах и библиотеках, да и собственноручно реализовать не составит особого труда. Между тем, опыт показывает, что, при всей своей простоте, метод начинает вызывать некоторые вопросы, когда возникает необходимость не просто посмотреть наличие дискреток в сигнале, но и выяснить их абсолютные значения, т.е. нормализовать полученный результат.

В этой статье я постараюсь объяснить, что же все-таки выдает в качестве результата fft (Fast Fourier transform) на примере MATLAB (и в качестве бонуса проведу небольшой ликбез по этому весьма полезному, на мой взгляд, языку).

MATLAB позволяет не заморачиваться с ручным удалением ненужных объектов, однако, при работе с более менее объемными массивами данных, имеет привычку капризничать и жаловаться на недостаток памяти. Для освобождения памяти используется процедура clear с указанием имени объекта, который необходимо удалить.

С этого и начнем. Так как все необходимое мы сгенерируем самостоятельно, можно смело удалять все, что накопилось в рабочем пространстве за активную сессию, просто добавив ключевое слово all:

clear all% Очистка памяти

Итак, прежде всего, зададим исходные данные для нашей модели. Фурье анализ идеально подходит для выделения гармонических сигналов на фоне помех. Для того чтобы продемонстрировать это, возьмем в качестве сигнала сумму некоторой постоянной и двух синусоид с разной частотой и амплитудой. Дисперсию шума возьмем в 3 раза больше амплитуды первой синусоиды. Так же зададим количество частотных полос, которые должен будет посчитать fft алгоритм. Точка с запятой в конце строк не обязательна, и если ее не ставить, результат вычисления функций и задания переменных будет дублироваться в командную строку, что можно использовать для отладки кода (однако, 512 значений сплошным полотном в командной строке вряд ли помогут вам, тем более что их вывод тоже занимает некоторое количество времени, так что все же лучше не забывать закрывать строки).

%% Параметры
Tm=5;% Длина сигнала (с)
Fd=512;% Частота дискретизации (Гц)
Ak=0.5;% Постоянная составляющая (Попугаев)
A1=1;% Амплитуда первой синусоиды (Попугаев)
A2=0.7;% Амплитуда второй синусоиды (Попугаев)
F1=13;% Частота первой синусоиды (Гц)
F2=42;% Частота второй синусоиды (Гц)
Phi1=0;% Начальная фаза первой синусоиды (Градусов)
Phi2=37;% Начальная фаза второй синусоиды (Градусов)
An=3*A1;% Дисперсия шума (Попугаев)
FftL=1024;% Количество линий Фурье спектра


MATLAB (Matrix Laboratory), как следует из названия, предназначен прежде всего для работы с массивами, практически все алгоритмы счета в нем оптимизированы для работы с векторами. Обилие удобных инструментов работы так же ненавязчиво подталкивает представлять как можно больше исходных данных в виде матриц. В частности, можно легко сгенерировать массив возрастающих (убывающих) величин с заданным шагом (1\Fd в данном примере):

%% Генерация рабочих массивов
T=0:1/Fd:Tm;% Массив отсчетов времени


Случайный Гауссов шум задается функцией randn, результатом которой является массив размерности, заданной в ее параметрах. Для единообразия зададим его в виде строки (первый параметр 1) длиной соответствующей длине нашего массива отсчетов времен, что в свою очередь вычислим функцией length.

Noise=An*randn(1,length(T));% Массив случайного шума длиной равной массиву времени

Символ * используется для обозначения перемножения. Т.к. чаще всего действия производятся над векторами, то и умножение подразумевается векторное, но так же легко можно, не перегружая этот оператор, использовать его для поэлементного перемножения, добавив перед ним точку (.*). При умножении вектора на скаляр точка перед умножением не является обязательной. Также добавленная точка сделает поэлементным возведение в степень и деление.

Signal=Ak+A1*sind((F1*360).*T+Phi1)+A2*sind((F2*360).*T+Phi2);% Массив сигнала (смесь 2х синусоид и постоянной составляющей)

Теперь перейдем к тому, ради чего и затевалась данная статья — функции fft(). Аргументами стандартной функции MATLAB являются сам сигнал (в нашем случае Signal), размерность вектора-результата (FftL), а также измерение.
Последний аргумент определяет вдоль какого измерения располагается сигнал в случае если на вход подается многомерный массив (Иногда этот параметр ошибочно принимают за размерность преобразования Фурье, но это не так. Хотя в MATLAB есть реализации 2-хмерного fft2() и многомерного fftn() алгоритмов). Так как наш сигнал представляет собой вектор, то его вполне можно опустить.
Рассмотрим сначала сигнал без примеси шума. В качестве результата мы получим вектор комплексных чисел. Это и есть представление нашего сигнала в частотном домене в показательной форме. Т.е. модули этих комплексных чисел представляют амплитуды соответствующих частот (точнее полосы частот см. дальше), а аргументы – их начальные фазы. И если полученная фаза, однозначно вычисляется в радианах, то с амплитудой и частотами не все так просто.
Например, если мы просто применим к нашему сигналу преобразование Фурье, и возьмем абсолютные значения вектора на выходе, то получим приблизительно следующую картинку:

image

Для построения двухмерных графиков удобно использовать функцию plot. Основные параметры, используемые в этой функции – одномерные массивы точек, первый задает ось ординат, второй – значение функции в соответствующих точках. Если передать только один массив, то он будет отображен с фиксированным шагом 1.
Если присмотреться к полученной картинке выяснится, что она несколько отличается от наших ожиданий. На приведенном графике 5 пиков вместо ожидаемых 3х (постоянная + 2 синусоиды), их амплитуды не совпадают с амплитудами исходных сигналов, и ось абсцисс вряд ли отображает частоты.

Прежде всего, следует учитывать, что счет алгоритма устроен таким образом, что перебираются не только положительные, но и отрицательные частоты и правая часть графика является «зеркальным» отображением реального спектра. Т.е. на самом деле 0 (которому соответствует постоянная часть сигнала) должен приходиться на середину массива. Ситуацию можно поправить, совершив циклический сдвиг на половину длины массива. Для этих целей, в MATLAB существует функция сдвига fftshift(), смещающая первый элемент в середину массива:

image

Теперь обратим внимание на ось значений.
Согласно теореме отсчетов (так же известной как теорема Найквиста-Шеннона или более патриотично теорема Котельникова) спектр дискретного сигнала будет ограничен половиной частоты дискретизации (Fd). Или в нашем случае –Fd/2 слева и Fd/2 справа. Т.е. весь полученный массив покрывает Fd частот. Отсюда, учитывая что мы знаем (вернее даже самостоятельно задали в качестве параметра) длину массива, получим частоты в виде массива значений от –Fd/2 до Fd/2 с шагом Fd/FftL (на самом деле крайняя правая частота будет меньше границы на один отсчет т.е. Fd/2-Fd/FftL):

image

Если посмотреть на фазы частот, можно заметить, что они равны отрицательным фазам соответствующих отрицательных частот. Учитывая равенство амплитуд левой и правой частей спектра и соответствие их фаз с точностью до знака, весь спектр будет эквивалентен своей положительной части с удвоенной амплитудой. Исключение составляет только 0 элемент, который не имеет зеркальной половины. Таким образом, можно избавиться от «непонятных» и зачастую ненужных отрицательных частот. Это можно было сделать сразу просто отбросив конец исходного массива и помножив оставшиеся элементы на 2 (за исключением постоянной составляющей):

image

Теперь это уже похоже на тот результат, который мы ожидаем. Единственное, что смущает теперь – это амплитуды. С этим все достаточно просто. Т.к. быстрое преобразование Фурье фактически представляет собой суммирование сигнала перемноженного на ядро преобразования (комплексную экспоненту) для каждой из частот, то реальный результат будет меньше полученного ровно в количество суммирований (частот в результате), т.е. полученный результат надо разделить на количество элементов в результате (не забываем, что под результатом понимается весь ответ, вместе с отброшенной частью, т.е. наше заданное FftL):

image

Стоит упомянуть еще одну вещь. В спектральном представлении вычисляется не значение сигнала на той частоте на которую попал алгоритм (как мы помним частоты следуют с шагом Fd/FftL), а значение в полосе (шириной равной шагу). Т.е. если в эту полосу попало несколько дискреток, то они суммируются. Для примера можно уменьшить количество линий в результате работы алгоритма:

image

Однако это не означает, что стоит сразу бездумно увеличивать точность работы, т.к. это тоже приводит к негативным последствиям, т.к. если разрешение будет сопоставимо с частотой дискретизации сигнала, в спектр полезут гармоники «окна», которые имеют отношение не к реальному сигналу, а к его дискретному представлению:

image

Или более близко окрестности одной из дискреток:

image

Код для нормировки fft будет выглядеть приблизительно следующим образом:

%% Спектральное представление сигнала
FftS=abs(fft(Signal,FftL));% Амплитуды преобразования Фурье сигнала
FftS=2*FftS./FftL;% Нормировка спектра по амплитуде
FftS(1)=FftS(1)/2;% Нормировка постоянной составляющей в спектре
FftSh=abs(fft(Signal+Noise,FftL));% Амплитуды преобразования Фурье смеси сигнал+шум
FftSh=2*FftSh./FftL;% Нормировка спектра по амплитуде
FftSh(1)=FftSh(1)/2;% Нормировка постоянной составляющей в спектре


Нам осталось только вывести результаты. Функция subplot позволяет разбить окно на несколько областей для отображения графиков.

%% Построение графиков
subplot(2,1,1);% Выбор области окна для построения
plot(T,Signal);% Построение сигнала
title('Сигнал');% Подпись графика
xlabel('Время (с)');% Подпись оси х графика
ylabel('Амплитуда (Попугаи)');% Подпись оси у графика
subplot(2,1,2);% Выбор области окна для построения
plot(T,Signal+Noise);% Построение смеси сигнал+шум
title('Сигнал+шум');% Подпись графика
xlabel('Время (с)');% Подпись оси х графика
ylabel('Амплитуда (Попугаи)');% Подпись оси у графика

F=0:Fd/FftL:Fd/2-1/FftL;% Массив частот вычисляемого спектра Фурье
figure% Создаем новое окно
subplot(2,1,1);% Выбор области окна для построения
plot(F,FftS(1:length(F)));% Построение спектра Фурье сигнала
title('Спектр сигнала');% Подпись графика
xlabel('Частота (Гц)');% Подпись оси х графика
ylabel('Амплитуда (Попугаи)');% Подпись оси у графика
subplot(2,1,2);% Выбор области окна для построения
plot(F,FftSh(1:length(F)));% Построение спектра Фурье сигнала
title('Спектр сигнала');% Подпись графика
xlabel('Частота (Гц)');% Подпись оси х графика
ylabel('Амплитуда (Попугаи)');% Подпись оси у графика


Результат выполнения кода будет выглядеть следующим образом:

image

image

Несмотря на то, что полезного сигнала не видно на фоне шума, спектральная характеристика позволяет определить его частоту и амплитуду.

Надеюсь данный текст был вам полезен.
Поделиться публикацией
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Реклама
Комментарии 52
  • +2
    выложи где-нибудь входящие данные
    • +4
      m-файл с полным скриптом:
      FScr.m
      • 0
        файл данных, с которыми работаете.
        • +1
          Сигнал генерирую в 19 строчке. В принципе, если хотите, могу записать синусоиду в WAV, но зачем?
          • 0
            то же самое в 1 команду


            • +1
              А еще pwelch, pmusic и т.д. Только топик был про нормировку вручную. MATLAB выбрал исключительно потому что без проблем визуализировать можно в одну команду.
              • 0
                Одно но: визуализировать лучше в логарифмическом масштабе (в децибелах). Куда как нагляднее получится. Как у коллеги двумя комментами выше.
                • +1
                  Выводил в линейном масштабе, чтобы явно было на входе 0.7 попугая и на спектре 0.7 попугая. А к какому основанию приводить попугаев в логарифмическом масштабе я не знаю.
                  • 0
                    Ну как вы верно заметили, измеряется все равно в попугаях, так что принципиальной разницы нет, важны соотношения между гармониками и уровень шума. Вообще в подобных случаях чаще всего используется единица dBfs — приведение к полной шкале. То есть уровень синусоидального сигнала амплитудой 1 (к примеру) берется за 0 dBfs, а все остальное — ниже.
                    • +1
                      Да к разному приводят… И к уровню тепловых шумов, и к единичным уровням, и к порогу слышимости… Да только моей целью было показать как БПФ выдает именно те попугаи, которые на входе.
      • +1
        Вы слабо разбираетесь в теме. Не пудрите голову людям.
        • +3
          Резко сказано… Но некоторые некорректности имеются.
          В качестве результата мы получим вектор комплексных чисел. Это и есть представление нашего сигнала в частотном домене в показательной форме. Т.е. модули этих комплексных чисел представляют амплитуды соответствующих частот (точнее полосы частот см. дальше), а аргументы – их начальные фазы.

          Во-первых, на графике показаны как раз модули этих комплексных чисел, во-вторых нет такого понятия «амплитуда полосы частот», в-третьих комплексное число не имеет аргументов — он сам есть число и может являться аргументов. Если имеются в виду действительная и мнимая части, то их две, а начальная фаза у синусоиды одна. Что имелось в виду?
          «следует учитывать, что счет алгоритма устроен таким образом, что перебираются не только положительные, но и отрицательные частоты»
          Отрицательных частот не существует, а симметрия графика — есть результат нормализации (возведения в квадрат).
          • +4
            Про график каюсь, надо бы указать, что строю именно амплитуды, собственно в самом скрипте видно, что я беру модуль.
            На счет амплитуды — не смог сформулировать по-другому, но имеется ввиду суммарный уровень частот попавших в полосу (при этом не могу сказать что там прямо таки сумма).
            Показательная форма представления комплексного числа — это не sin(x)-i*cos(x), а a*exp(i*phi), где а — модуль, а phi — аргумент, он же и будет начальной фазой (может стоит внести это пояснение в текст?).
            На счет «отрицательных» частот, тоже вопрос для дискуссии. Я БПФ представляю как последовательный перебор корреляции сигнала с синусоидами различных частот начиная с отрицательных индексов. Физического смысла оно действительно не имеет, но в рамках алгоритма вполне себе живет.
            • 0
              Я о том и говорю. Приветственная фраза обязывает излагать так, чтобы не вводить в заблуждение «юного читателя».
              • 0
                «Куда ставить-то?», в смысле, что поправить-то на ваш взгляд надо?
                • 0
                  Если кроме меня к формулировкам так дотошно никто не относится, то особо нечего менять. Тут больше вопрос стиля, наверное.
                  Как видно, несмотря на то, что сигнал весьма сильно скрыт шумом, в спектральном представлении все еще отчетливо видны дискретные составляющие (хоть и немного искаженные по амплитуде шумом).

                  Фраза «весьма сильно скрыт шумом» не очень хороша для технической публикации. Поэтично что ли… «Видны дискретные составляющие» — это тоже нехорошо. У вас же сам сигнал в дискретном времени, в нем каждый отсчет — дискретная составляющая. Легко ввести читателя в заблуждение, а знающим труднее понять. Можно же было так написать:
                  «Несмотря на то, что полезного сигнала не видно на фоне шума, спектральная характеристика позволяет определить его частоту и амплитуду.» Хотя и эта формулировка не безгрешна.
              • 0
                имеется ввиду суммарный уровень частот попавших в полосу (при этом не могу сказать что там прямо таки сумма

                Это обобщенно и проще можно назвать энергией в данной элементарной полосе частот.
                http://ru.wikipedia.org/wiki/%D0%90%D0%BD%D0%B0%D0%BB%D0%B8%D0%B7%D0%B0%D1%82%D0%BE%D1%80_%D1%81%D0%BF%D0%B5%D0%BA%D1%82%D1%80%D0%B0
                • 0
                  Да что такое… неудобно ссылки вставлять. Простите.
                  Анализатор спектра
                  • 0
                    Только тут спектр амплитудный, а не энергический (для наглядности картинок делал), поэтому про энергию в полосе частот тоже говорить не корректно на мой взгляд.
                    • 0
                      А собственно почему не энергетический. Как раз на мой взгляд нечто похожее на гистограмму (по сути, не по внешнему виду). Нужно еще заметить, что по этим графикам нельзя определить был ли полезный сигнал на протяжении всего анализированного интервала времени одинаков по амплитуде. Его амплитуда могла быть как постоянной, так и меняться во времени. Это же графики квадратичных величин. Как раз энергии и есть. Но этим графикам исходный сигнал не восстановишь в первозданном виде. Каждый пик по сути это энергия соответствующего интервала частот (гармоник), приведенная к общей энергии спектрограммы (сигнала). У меня в свое время поэтому руки не доходили ответа на вопрос «в каких попугаях» построены эти графики. Потом вопрос утратил актуальность.
                      • 0
                        Это именно амплитудный спектр, если вы обратите внимание амплитуды пиков в спектре совпадают с заданными в качестве параметров. А если обратить внимание еще и на фазовый спектр (его не приводил в связи с непрезентабельностью), то можно в точности восстановить весь сигнал на протяжении всех 5 секунд. Собственно это и хотел показать, что пики находятся именно на своих местах и принимают реальные значения сигнала во временном домене безотносительно к типу сигнала (паскали, вольты, метры). А на счет изменения амплитуды со временем, так это уже модуляция, а спектр модулированного сигнала будет отличаться наличием гармоник модулирующего и это уже будет предметом анализа о чем речи в топике не идет.
                        • +1
                          А можно линк на тему того, как огибающая (кривая амплит. модуляции) появляется на спектрограмме?
            • +1
              Отрицательные частоты имеют вполне такой себе физический смысл. И в рамках очень простого эксперимента можно реально увидеть их на спектро-анализаторе. Возьмите и подайте на спектроанализатор последовательность прямоугольных видеоимпульсов, на картинке вы увидите спектр этой последовательности который будет представлять из себе 8ю картинку из этой статьи с максимальной частотой равной 0, и соответственно видеть мы будем только половину этого графика. Если мы начнем заполнять наш видеоимпульс несущей частотой то картинка, с ростом частоты начнет смещаться вправо по оси х. При несущих частотах близких к нулю левая от максимума(несущей частоты) часть графика будет отличаться от правой. Это потому, что в этой области частот происходит наложения двух частей спектра сигнала — спектра с положительной центральной частотой и спектра с отрицательной центральной частотой, правый хвост последнего (боковые лепестки) вылазит в область спектра положительных частот и в результате наложения левой части спектра с положительной центральной частотой и правой части спектра с отрицательной несущей частотой мы увидим не симметричный относительна несущей частоты спектр, а искаженный.
              • 0
                Вот как я и говорил — ввели в заблуждение…
                Задача для Буратино:
                Малыш, представь, что у тебя в кармане три яблока и тут некто взял у тебя одно яблоко. Сколько у тебя в кармане яблок? Не торопись с ответом и учти, что в нашей школе принято отнимать от ответа «СТОПИЦЦОТ» штук.
                • 0
                  На всякий случай оговорюсь… Буратино — это я фильм детский вспомнил. Никого обзывать не собирался.
                  • 0
                    С Буратиной все ясно, но вот кого ввели в заблуждение я честно говоря не понял. Если меня то хотелось бы понять где я заблуждаюсь, да и в любом случае любые заблуждения нужно искоренять :)
                    • +1
                      Я про то, что линии на графике, строго говоря, не есть частоты, а лишь их изображение, интерпретация, визуализация. Какой же это физический смысл? То, что смещена ось частот на индикаторе не говорит об отрицательности частот. И если брать несущую частоту — ее близость к нулю (абсолютному, если не оговорено опять же смещение оси частот) означает близость несущего сигнала к постоянному. «Постояннее», чем постоянный сигнал, никакой сигнал быть не может. — Какой каламбур получился… Ужас.

                      На мой взгляд нельзя подменять парадигмы искусственно выведенными (производными) формулировками.
                      • 0
                        А я не говорил про графики. Я говорил про реальный спектроанализатор, и рассказал про эксперимент который позволяет ощутить реальность несуществующих отрицательных частот. Во как «реальность несуществующих»:)
                        • +1
                          «Наш преподаватель по физике, рассказывавший про эл/маг. волны, был настолько пьян, что мы реально увидели эти волны.
                          Какая разница какой индикатор? Растр/Вектор на экране монитора, подсвеченный люминесцентный экранчик осциллографа или любой другой?
                          • –1
                            Наверное тут имелась ввиду амплитудная модуляция, а не аддитивная смесь. Там реально возможна картинка с пиками, пляшущими вокруг несущей (и слева и справа). Только все же отрицательные частоты действительно не несут особого смысла, это тот же самый синус, но повернутый на 180 градусов.
                            • +1
                              Можете в википедии почитать про Отрицательные частоты . Если вращается в одну сторону, то имеет положительную частоту, а если в другую, то отрицательную.
                              Получается, что частота вдруг превратилась в векторную величину и была частота сама по себе, а стала некоей относительной или зависимой.
                              Как профессиональный жаргон сойдет.
                              • 0
                                Мы вообще-то говорим про сигнал. Он «вращается» только в одну сторону — вперед во времени, а вот колебаться с одной и той же частотой он может в разных фазах, например в противофазе:
                                sin(a*x+180)=-sin(a*x)=sin(-a*x)
                                Вот вам и формальная отрицательная частота.
                  • 0
                    «максимальной частотой равной 0» следует читать как «с максимумом в частоте равной нулю»
                    • 0
                      Вот об этом я и говорю. Нужно быть осторожнее в таких вопросах. Прочитайте снова эти две формулировки — одно и то же. Вы имели в виду «частоту, соответствующую максимальному значению графика»?
                      • 0
                        Да, из-за описки смысл текста мог поменяться. тут я накосячил.
                        • 0
                          Во-во… «Требую долива отстоя после пены» или «Казнить (,) нелья (!,) помиловать»…
                    • –1
                      > Отрицательные частоты имеют вполне такой себе физический смысл. И в рамках очень простого эксперимента можно реально увидеть их на спектро-анализаторе.

                      Круто! На спектро-анализаторе!

                      А в природе нельзя?
                    • 0
                      Не стал бы называть это некорректностями. Человек не знает что такое окно и зачем оно нужно в спектральном анализе. Код я смотрел нет там никакого окна. Кроме того, автор не понимает что спектр расплылся (при увеличении длины БПФ) не от «гармоник окна» а от некореектной (не кратной) передискретизации.

                      О таких мелочах как отсутствия связи с физическим временем, а как следствие некоррекнтое использование размерности (секунды, герцы) в изначально цифровом синале я молчу.

                      Всем минусующим карму пламенный привет.
                      • 0
                        Вы точно читали статью и смотрели код? Целью топика было не объяснение особенностей методов уменьшения влияния краевых эффектов, а попытка показать, как можно убедиться, что результатом выполнения кода является именно то, что вы ожидаете увидеть. Сигнал генерировался относительно смоделированного времени. С заданными частотой дискретизации и длиной окна (фактически 5 секундное прямоугольное окно, или быть может вы предполагали, что я сгенерировал бесконечный сигнал?). И частоты приводились именно через частоту дискретизации сигнала.

                        А про карму, рискну предположить что кому-то не очень понравился «конструктивный» подход к написанию комментариев. Со своей стороны я бы тоже попросил больше конкретики, если я где-то заблуждаюсь укажите конкретно где, думаю не мне одному будет интересно выяснить какие-то детали которые я не понимаю или упускаю.
                  • +1
                    Уххх как знакомо, у меня тема диплома была «Speech processing using DSP» — ну как то так, дословно уже не помню. Помучался я с этим ( pastebin.com/f12259e0d ) кодом…

                    x -> fft -> ifft -> profit
                    • 0
                      Ага. А еще забыл про фазы и на выходе кваказябры.
                      • 0
                        Я уже всё забыл :D
                    • +1
                      Добавьте пожалуйста в теги Фурье или преобразование Фурье, а то непонятно как люди искать потом этот топик будут, ведь главное тут не Матлаб, а алгоритм, о чём и говорит нам название блога =)
                    • 0
                      Где вы были 4 дня назад… Когда я лабу писал
                      • 0
                        Хорошая статья, спасибо. Давно собирался поковырять MATLAB, вот будет повод.
                        Было время, FFT на C++ делал для дипломной. Тоже, помню, вылезли эти зеркальные отображения, физический смысл их так и не понял.
                        • 0
                          А это и не физический смысл, а «побочный эффект» быстрого преобразования Фурье.
                        • 0
                          Попробовал распознать ноту Ми (82.406 Гц). Получилось.


                          clear all% Очистка памяти
                          [Signal,hc,bits] = wavread('E2.wav');%Файл с длинным звучанием ноты
                          
                          %% Параметры
                          Tm=0.5;% Длина сигнала (с)
                          Fd=hc;% Частота дискретизации (Гц), у меня была 22050 из wav характеритстика
                          FftL=Fd*2;% Количество линий Фурье спектра
                          Epsilon=50;% Окрестность максимальной частоты для наглядного графика n.3
                          Signal = transpose(Signal(1:(Fd*Tm)));% все же берем отрезок сигнала
                          
                          %% Спектральное представление сигнала
                          FftS=abs(fft(Signal,FftL));% Амплитуды преобразования Фурье сигнала
                          FftS=2*FftS./FftL;% Нормировка спектра по амплитуде
                          FftS(1)=FftS(1)/2;% Нормировка постоянной составляющей в спектре
                          
                          F=0:Fd/FftL:Fd/2-1/FftL;% Массив частот вычисляемого спектра Фурье
                          [C,i] = max(FftS);%Максимум и его аргумент вмассиве частот спектра
                          
                          %% Построение графиков
                          figure( 'NumberTitle','off',...
                                  'Name','BorisPlus: распознование ноты (Ми большой октавы/E2) из файла')
                          
                          subplot(3,1,1);% Выбор области окна для построения
                          plot(Signal);% Построение сигнала
                          title('Сигнал');% Подпись графика
                          xlabel('Время(с)');% Подпись оси х графика
                          ylabel('Ампл.');% Подпись оси у графика
                          
                          subplot(3,1,2);% Выбор области окна для построения
                          plot(F,FftS(1:length(F)));% Построение спектра Фурье сигнала
                          title('Спектр сигнала');% Подпись графика
                          xlabel('Частота(Гц)');% Подпись оси х графика
                          ylabel('Ампл.');% Подпись оси у графика
                          
                          subplot(3,1,3);% Выбор области окна для построения
                          plot(   F(max(1,i-Epsilon):min(length(F),i+Epsilon)),...
                                  FftS(max(1,i-Epsilon):min(length(F),i+Epsilon)),'-r',...
                                  F(i),FftS(i),'ko');% Посмотрим наглянее у максимума
                          title('Спектр сигнала в окрестности максимума');
                          xlabel('Частота(Гц)');% Подпись оси х графика
                          ylabel('Ампл.');% Подпись оси у графика
                          text(F(i)+1,FftS(i),[num2str(F(i)) ' ' 'Гц']);% Подпись точки максимума
                          

                          Есть мысль развить сие. Спасибо за толчок.
                          • 0
                            Рад что кому-то полезно.
                          • 0
                            Ну а теперь, чтоб было все просто красиво, как с подбором длинны иглы в опыте Бюффона, чтоб получить точное зачение Пи за одно-два бросания. Распознаем ноту Ля большой октавы (110 Гц). О-па! Чётко.
                            • 0
                              Лайк почему-то прожать не могу, но спасибо тебе, друг!

                              Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.