Еще один способ разложения сигнала в спектр

Привет всем, здесь я хочу рассказать про алгоритм анализа звукового сигнала, позволяющий разобрать сигнал на отдельные волны, конечно 100% точности он не дает, но тем не менее результат на мой взгляд довольно неплох.

image

Лучше всего видна работа на какой-нибудь музыке:


И ссылки на другие примеры разных жанров. Metaldeth-Tornado of souls:


Out of Touch:


Итак для разложения нужно сделать следующие шаги:

— Из исходного сигнала нужно получить 8 промежуточных сигналов;
— Из этих промежуточных сигналов и исходного сигнала нужно получить 8 сигналов — слоев, которые можно будет разобрать на отдельные волны;
— Посчитать сколько в каждом слое волн и какая у них амплитуда.

Теперь подробнее о каждом этапе: для того чтобы получить промежуточный сигнал нужно взять производную от исходного сигнала. По сути это производная дискретной функции. Чтобы найти ее для каждого момента исходного сигнала нужно задать 1 параметр: период за который эта производная находится. Значение производной это коэффициент наклона в заданном интервале, можно найти например методом наименьших квадратов.

Требуется вычислить 8 промежуточных сигналов с 8 разными периодами. Самый простой набор периодов: 4, 8, 16, 32, 64, 128, 256, 512. Когда период задан, для каждого отсчета сигнала вычисляется производная по формуле наименьших квадратов. Это как скользящее среднее, только здесь не скользящее среднее а скользящая производная текущего интервала.

Таким образом получается 8 производных сигналов и 1 исходный. Теперь каждый из 8 производных сигналов нужно проинтегрировать. В данном случае это значит что каждый следующий семпл равен сумме всех предыдущих семплов. После этого получается 8 промежуточных слоев.

Следующий шаг — получение слоев, которые можно будет разобрать на отдельные волны. Итак теперь нужно получить 8 слоев. Слои вычисляются так:

слой0=промежуточный0-исх сигнал
слой1=промежуточный1-промежуточный0
слой2=промежуточный2-промежуточный1
слой3=промежуточный3-промежуточный2
слой4=промежуточный4-промежуточный3
слой5=промежуточный5-промежуточный4
слой6=промежуточный6-промежуточный5
слой7=промежуточный7-промежуточный6
слой8=промежуточный7

Последний слой является не разницей а просто равен последнему промежуточному сигналу.

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

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

image

Код который вычисляет промежуточный сигнал выглядит вот так:

здесь wavesize – число семплов
signal[] – массив с исходным сигналом
SY=0,SX=0,SXX=0,SXY=0,Ky=0 — переменные типа float
Step2=STEP/2 где STEP это период (4,8,16,32,64,128,256,512)

for(int i=Step2;i<wavesize-Step2;i++){
	SY=0,SX=0,SXX=0,SXY=0,Ky=0;

	for(int j=i-Step2,fromZ=0;j<i+Step2;j++,fromZ++){
		SX+=fromZ;
		SY+=signal[j];
		SXX+=fromZ*fromZ;
		SXY+=fromZ*signal[j];
	}
	Ky=float((STEP)*SXY-SX*SY)/float((STEP)*SXX-SX*SX);
		
	OutSignal[i]=OutSignal[i-1]+Ky;	
}

Чтобы вычесть один сигнал из другого достаточно просто вычесть каждый семпл из каждого.
Например для 0 слоя:

for(int i=0;i<wavesize-1;i++)		
 layer0[i]=OutSignal0[i]-Signal[i];

Если сложить все слои, и взять последний с обратным знаком то получится исходный сигнал, таким образом умножая какой-нибудь слой, возможно сделать фильтр частот. Далее вопрос заключается в том как считать амплитуды конкретных гармоник. Дело в том что в постоянном интервале, например = 4000 семплов, может быть очень много коротких волн и сравнительно мало длинных.

Можно конечно найти средние амплитуды для каждого слоя и сложить. Но этот способ не очень хорош так как длинных волн очень мало, а их амплитуда обычно очень большая, и получается сильная неравномерность амплитуды в сторону низких частот.

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

Вообще, не думаю что человек воспринимает звук как разложение в спектр, скорее звук состоит из звуковых образов, которые состоят из волн разной длинны. Соответственно громкость какого либо звука это скорее средняя громкость всех волн из которых звук состоит. Но пока что не понятно какие параметры составляют звуковой образ, возможно средняя частота, стандартное отклонение или что то еще.
Поделиться публикацией
Ой, у вас баннер убежал!

Ну. И что?
Реклама
Комментарии 42
  • 0
    Вы не пробовали 3-4 заранее известных синусоиды смиксовать и подать и посмотреть что получится? Очень интересно, если можете, добавьте в статью.
  • 0
    То что вы описали очень похоже на вейвлет Хаара, нет? Каждый следующий шаг раскладывает сигнал на низкочастотную и высокочастотную составляющие.
    • 0
      А что именно похоже? Посмотрел бегло про преобразование Хаара, оно вообще не про спектральный анализ. Теперь посмотрел подробнее, действительно есть упоминания про частотный фильтр, но там насколько я понял другая идея
      • 0
        Хаар вообще-то математик и он много чего понапридумывал. При чем тут преобразование? Я ж вам дал ссылку на вейвлет — очень похоже на ваш способ.
        • 0
          Посмотрел по ссылке пункт: преобразование для одномерного сигнала. Там исходный сигнал раскладывается на 2 других, но чтобы их получить нужно складывать/вычитать соседние точки. Речи про производные там нет.
          • +1
            Да, мне тоже вспомнилось дискретное вейвлет-преобразование Хаара. Спектр сигнала при этом вычисляется изменением ширины вейвлет-функции: 2, 4, 8 и т.д. Функция умножается на исходный сигнал, получается «послойное» разложение сигнала на частоты. Может быть не в чистом виде так, но очень похоже.
            • 0
              По вашей ссылке есть код. Сравните с тем что в статье.
              • +2
                Код — всего лишь способ реализации идеи (алгоритма), одну идею можно реализовать разным кодом, нет смысла сравнивать коды напрямую.
                Я же написал, что ваша идея мне напоминает вейвлет-преобразование и описал какое именно.
                Я не утверждал, что у вас вейвлеты, всего лишь сказал, что ваша идея похожа на вейвлеты Хаара.
                Так это или нет решать вам, вы лучше меня представляете вашу же идею выделения спектра и её реализацию в виде программы.
    • +4
      Зачем вся эта цветомузыка? Вместо неё должен быть пример спектрограммы, полученной с помощью алгоритма. Можно в сравнении с классическим STFT.
      • +2
        Мне показалось, или Вы пере-изобрели Преобразование Фурье?
        • –1

          Автор по крайней мере попытался его переизобрести, но куда-то потерял фазу.
          Действительно надо посмотреть на результат работы с одной/несколькими синусоидами.

          • –1
            Почему Вы так думаете?
          • –1
            Здесь же Гиктаймс а не хабр, а эта цветомузыка и есть спектрограмма.(это ответ для Sadler)
            • +2
              Претензия к Вам лишь в том, что вы выбрали крайне неинформативное представление результатов, и я не знаю, продиктовано это отсутствием опыта или желанием скрыть общую картину. На хабре Вы бы ещё добавили код этого безобразия.
            • 0
              Вообще, не думаю что человек воспринимает звук как разложение в спектр, скорее звук состоит из звуковых образов, которые состоят из волн разной длинны. Соответственно громкость какого либо звука это скорее средняя громкость всех волн из которых звук состоит. Но пока что не понятно какие параметры составляют звуковой образ, возможно средняя частота, стандартное отклонение или что то еще.

              Человек воспринимает звук субъективно. То есть по сути мы понятия не имеем как конкретный человек слышит тоже само что и мы. Описывает он это теми же самыми словами, но что по факту, мы можем только гадать.
              Громкость в ухе передается тактом от 20 до 100 гц.
              А про звуковые образы и обобщение различных сенсоров (ака сенсорные нейроны, которых миллионы), можно почитать кучу книг, и написать еще не меньшую кучу. :)
              • 0
                Да, субьективно. Но все различают скачки (удары) и монотонный звук. И также различают гласные, согласные. Т.е мозг таки имеет(учится иметь) стандартные критерии в сигнале, на основании которых делает вывод что это за звук.
                А про передачу громкости в ухе, имеется ввиду частота спайков нейронов?
              • 0
                1. Этап до получения промежуточный0… промежуточный8 очень долго описывается, а на самом деле это просто пачка ФНЧ. Если это понять — и проследить алгоритм (предсказав его поведение) проще, и ускорить можно.
                2. Как вам посоветовали, попробуйте на смеси двух синусоид. Сразу скажу, тут будет косяк: если обе частоты попадают в один слой — ваш детектор выдаст одну частоту (с переменной амплитудой).
                3. Хотелось бы видеть сравнение по быстродействию с БПФ.

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


                Бонус (раз уж цветомузыка): фильтр типа y(t) = 0.9y(t-T)+0.1x(t) выделяет частоту 1/T и её гармоники, если сделать пачку таких фильтров для 12 нот какой-то октавы — можно отрисовать 12 красивых осциллограмм, дёргающихся в такт нотам.

                • 0
                  1)Ну конечно это пачка фильтров. Но ведь все дело именно в том как фильтр работает.

                  2)Про смесь 2 синусоид, если сложить 2 синусоиды частоты которых отличаются меньше чем в 2 раза то никакой алгоритм не даст 2 этих синусоиды. Поправьте если не так.

                  По поводу фильтра типа y(t) = 0.9y(t-T)+0.1x(t), а есть примеры работы этого фильтра которые можно посмотреть?
                  • 0
                    2 — преобразование Фурье даст (оно для этого и предназначено)
                    Пример — с ходу не найду, если покопаться в _очень старых завалах_ — возможно, накопаю исходник на паскале лохматого года под DOS, цифрующий звук с Sound Blaster Pro и отображающий 12 графиков в 13h видеорежиме VGA (320*240, 256 цветов :-))
                    Наверное, это уже даже не запустить без DosBox… Но визуально эффект красивый, нота в мелодии отображается сигналом на одном из графиков, причём форма сигнала — не синусоида (есть гармоники)

                    Да, регулировка «чувствительности» (добротности фильтров) — баланс между 0.9 и 0.1. Сейчас уже не вспомню — возможно, у меня было 0.95 и 0.05, или как-то так.
                    • 0
                      2 — но только я сомневаюсь что оно даст 2 чистые синусоиды, скорее весь интервал между частотами будет не 0 (я пробовал юзать FFT раньше, так что примерно представляю какой оно выдает результат).

                      А как называется этот фильтр?
                      • 0
                        Чистые синусоиды (точнее, амплитуды/фазы двух синусоид) — даст. Единственное — если окно не кратно частоте синусоиды, она ещё на какие-то частоты залезает (см. «боковые лепестки», «растекание спектра»), лечится это с помощью оконной функции.

                        Фильтр — гребенчатый фильтр или что-то типа того. Крайне примитивная и быстро считающаяся вещь, но именно для красивого отображения нот хорош.
                        • 0
                          Однозначно, различит. Предел различимости частот 1/T, где Т — длительность сэмпла. В два раза выше должна быть частота взятия отсчётов, чем высшая измеряемая частота в спектре. А практически, для качественных измерений частота взятия отсчётов берётся в 3 раза выше. Я как практик говорю, сам разрабатываю.
                          • 0
                            Да тут речь не про теорему Найквиста/Котельникова, а про результат после анализа. Т.е на мой взгляд после анализа суммы двух близких частот F1/F2<2, FFT даст весь набор между F1 и F2, возможно у F1 и F2 будут больше амплитуды, но будет и много лишнего.
                            • 0
                              нет, результат довольно чистый. Обязательно накладывается оконная функция. чтобы, как написали выше не было «растекания спектра». И результат настолько чистый, что побочные гармоники не наблюдаются практически.
                              • 0
                                Если взять например, волну с одной единственной частотой, но периодически меняющейся амплитудой, то FFT дает вот что:image
                                вот как выглядит сигнал: image
                                Если конечно windows player корректно юзает FFT
                                • 0
                                  Странно и нелепо. Вы бы написали простенькую программу да и анализированли. вот FFT отработало на массиве 16к точек, три синусоиды с частотами 1000, 1004, 1008

                                  Даже окно не пришлось накладывать, но для реальных сигналов лучше наложить, использую окно Хэмминга если интересно.
                                  • 0
                                    Ведь этот алгоритм написан уже много раз, зачем писать еще 1? Можете дать ссылку на программу которая на ваш взгляд корректно отрабатывает FFT. Например вот эта программа норм: ua3vvm.qrz.ru/gram.htm?
                                    • 0
                                      Можете дать ссылку на программу которая на ваш взгляд корректно отрабатывает FFT.

                                      Я могу вам дать программу
                                      // Pow - степень массива сэмплов, т.е. для 16к Pow = 14. 
                                      // Результат будет симметричен в этом случае относительно 8к, т.е. половину массива сэмпла
                                      void DP__Complex_FFT(double * Samples_Source,
                                                           double * Desination_R,
                                                           double * Desination_I,
                                                           u4x Pow)
                                      {
                                      
                                      	int Le, Le1, jt, it, ip, lt, kt, num;
                                      	double re_T, re_U, re_W, im_T, im_U, im_W;
                                      
                                      	num = 1 << Pow;
                                      	for (it = 0; it < num; it++) {Desination_R[it] = Samples_Source[it]; Desination_I[it] = 0;}
                                      
                                      	for (it = 0; it < num; it++)
                                      	{
                                      		kt = num / 2;
                                      		jt = 0;
                                      		Le = 1;
                                      		Le1 = it;
                                      		for (lt = 1; lt <= Pow; lt++)
                                      		{
                                      			if (kt <= Le1)
                                      			{
                                      				jt = jt + Le;
                                      				Le1 = Le1 - kt;
                                      			}
                                      			Le = Le * 2;
                                      			kt = kt / 2;
                                      		}
                                      
                                      		if (it < jt)
                                      		{
                                      			re_T = Samples_Source[jt];
                                                              im_T = 0;
                                      			Desination_R[jt] = Desination_R[it];
                                                              Desination_I[jt] = Desination_I[it];
                                      			Desination_R[it] = re_T;
                                                              Desination_I[it] = im_T;
                                      		}
                                      	}
                                      
                                      	for (lt = 1; lt <= Pow; lt++)
                                      	{
                                      		Le = 1;
                                      		for (jt = 1; jt <= lt; jt++) Le = Le * 2;
                                      		Le1 = Le / 2;
                                      		re_U = 1.0;
                                      		im_U = 0.0;
                                      		re_W = cos(M_PI / Le1);
                                      		im_W = sin(M_PI / Le1);
                                      		for (jt = 0; jt < Le1; jt++)
                                      		{
                                      			it = jt;
                                      			while (it < num - Le1)
                                      			{
                                      				ip = it + Le1;
                                      				re_T = Desination_R[ip] * re_U - Desination_I[ip] * im_U;
                                      				im_T = Desination_R[ip] * im_U + Desination_I[ip] * re_U;
                                      				Desination_R[ip] = Desination_R[it] - re_T;
                                      				Desination_I[ip] = Desination_I[it] - im_T;
                                      				Desination_R[it] = Desination_R[it] + re_T;
                                      				Desination_I[it] = Desination_I[it] + im_T;
                                      				it = it + Le;
                                      			}
                                      			re_T = re_U * re_W - im_U * im_W;
                                      			im_T = re_U * im_W + im_U * re_W;
                                      			re_U = re_T;
                                      			im_U = im_T;
                                      		}
                                                      
                                      	}
                                      
                                      }
                                      
                                      
                                      
                                      
                                      
                                      
                                      
                                      // Модуль комплексного спектра
                                      // Рабочие массивы I и R - превращаются в Desination
                                      
                                      double DP__Amplitude_FFT(double * Source_R,
                                                              double * Source_I,
                                                              double * Desination,
                                                              u4x Steps_amount)
                                      {
                                      
                                        double Max = 0.0;
                                      
                                        double Y;
                                      
                                        for(u4x i = 0; i<Steps_amount; i++)
                                        {
                                          
                                          Desination[i] = sqrt(((double)Source_R[i])*((double)Source_R[i]) + ((double)Source_I[i])*((double)Source_I[i]));
                                      
                                          if(Desination[i] > Max)
                                      
                                            Max = Desination[i];
                                          
                                        }
                                      
                                        return(Max);
                                      
                                      }


                                      И вообще это иллюстрация к тому, что при FFT не вылазит всяких артефактов, которых вы опасаетесь.
                                      Т.е на мой взгляд после анализа суммы двух близких частот F1/F2<2, FFT даст весь набор между F1 и F2, возможно у F1 и F2 будут больше амплитуды, но будет и много лишнего.
                                      • 0
                                        Ок, спасибо, на самом деле результат на Вашем скриншоте выглядит круто. Повожусь с кодом который Вы дали.
                                        • 0
                                          int main(int argc, char *argv[])
                                          {
                                              QCoreApplication a(argc, argv);
                                              printf("Start analyze\n");
                                              uint size=16384;
                                              double* Signal=new double[size];
                                              double* Desination_R=new double[size];
                                              double* Desination_I=new double[size];
                                              double* Desination=new double[size];
                                              memset(Signal,0,sizeof(double)*size);
                                              memset(Desination_R,0,sizeof(double)*size);
                                              memset(Desination_I,0,sizeof(double)*size);
                                              memset(Desination,0,sizeof(double)*size);
                                          
                                          
                                              for(uint i=0;i<size;i++)
                                                  Signal[i]=1000*sin(2*3.1415*float(i)*1000.0/44100.0)+1000*sin(2*3.1415*float(i)*1004.0/44100.0)+1000*sin(2*3.1415*float(i)*1008.0/44100.0);
                                          
                                          
                                              DP__Complex_FFT(Signal,Desination_R,Desination_I,14);
                                              //DP__Amplitude_FFT(Desination_R,Desination_I,Desination,14);
                                          
                                              for(uint i=0;i<size;i++){
                                                 // if(Desination_R[i]>1)
                                                 //     printf("%i  %f\n",i,Desination_R[i]);
                                                  if(Desination_I[i]>1)
                                                      printf("%i  %f\n",i,Desination_I[i]);
                                          
                                              }
                                          
                                          
                                              return a.exec();
                                          }


                                          Что то получается какая-то фигня у меня. Что в массиве R что в массиве I
                                          • 0
                                            нужно получать DP__Amplitude_FFT, и там уже спектр будет
                                            • 0
                                              пробовал и с DP__Amplitude_FFT (если ее раскоментировать и посмотреть содержимое Desination, то получаются какие то огромные числа)

                                              for(uint i=0;i<size;i++)
                                                      Signal[i]=1000.0*sin(2.0*3.1415*float(i)*1000.0/44100.0)+1000.0*sin(2.0*3.1415*float(i)*1004.0/44100.0)+1000.0*sin(2.0*3.1415*float(i)*1008.0/44100.0);
                                              
                                              
                                                  DP__Complex_FFT(Signal,Desination_R,Desination_I,14);
                                                  DP__Amplitude_FFT(Desination_R,Desination_I,Desination,size);
                                              
                                                  for(uint i=0;i<size;i++){
                                                      if(Desination[i]>1)
                                                         printf("%i  %f\n",i,Desination[i]);


                                              • 0
                                                может я неправильно обозначил u4x, я так понял что это целый тип, обозначил его как short. #define u4x short
                                                • 0
                                                  Большие числа у вас потому что большие амплитуды синусов, оставьте 1.
                                                  • 0
                                                    Поставил амплитуду 1, в результате получилось куча значений больше 1, даже есть пара значений больше 1000. Правильно ли я понимаю что в массиве Desination[i], амплитуда гармоники — значение Desination[i], а частота — индекс i?
                                                    • 0
                                                      Попробуйте одну синусоиду,
                                                      Signal[i]= sin(2.0*3.1415*float(i)*1000.0/16384);

                                                      и почему вы делаете так
                                                      if(Desination[i]>1)?

                                                      про Desination[i] понимаете верно
                                                      • 0
                                                        Все, заработало. Потестил разные варианты, работает очень даже норм.
                                                        Desination[i]>1, это чтобы не выводить 0 значения.
                            • 0
                              s/240/200/
                              (на случай, если придёт кто-то, кто помнит историю, и заметит мою ошибку)
                        • 0
                          Скажите, а вы не пробовали замерять скорость работы вашего алгоритма и быстрых преобразований Фурье? А также сравнить результаты работы вашего алгоритма и БПФ?
                          • 0
                            Не пробовал, да и та реализация что есть сейчас, далеко не самая быстрая.
                          • 0
                            И что касается фильтров и что касается преобразований Фурье, все они анализируют определенный временной участок. Чем более крутой фильтр, выше его порядок, тем дольше он «устанавливается». То-же относится и к Фурье. Преобразование Фурье дает среднеее значение амплитуды и фазы анализируемого фрагмента. Чем больший участок времени вы берете, тем точнее разделяются частоты, но за все надо платить, и вычисленные вариации амплитуды сглаживаются. Вообще, преобразование Фурье применяется для анализа бесконечной волновой функции, речь и музыка к таковым не относятся. Выше кто-то исследовал работу Фурье при изменении амплитуды. Господа, вы забываете, что изменение амплитуды, т.е амплитудная модуляция, приводит к изменению амплитуд и соседних частот.
                            Если вы берете какой-то фрагмент для анализа, то в него попадает один период самой низкой частоты и N/2 периодов самой высокой частоты. Естественно, вариативность нижних и верхних частот различаются. Хотите на каждом промежутке времени знать максимально мгновенное значение амплитуды гармоники — считайте Фурье за один ее период. И мир не ограничен только БФП. Есть и классический вариант
                            SinA = 1/N (Sum (X(i)*Sin (2*Pi*k*i/N))
                            CosA = 1/N (Sum (X(i)*Cos (2*Pi*k*i/N))
                            Ak = Sqr(SinA^2+CosA^2), где k — номер гармоники, N размер массива, Sum считается по массиву X с индексом i от 0 до N-1
                            И не забывайте накладывать окна (Хемминга например). Преобразования Фурье очень чувствительны к выбросам на концах анализируемого массива.

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

                            Самое читаемое