Лучше всего видна работа на какой-нибудь музыке:
И ссылки на другие примеры разных жанров. 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 вариант.
Теперь чтобы разобрать слои на отдельные волны нужно просто посчитать участки где значения увеличиваются и где они уменьшаются. Собственно длительность участков и есть их длины волн. Участки сигнала где значения сигнала постоянны нужно просто пропустить. Чтобы найти амплитуду сигнала в спектре в некотором интервале, нужно сложить все амплитуды волн умноженные на их длину.
Код который вычисляет промежуточный сигнал выглядит вот так:
здесь 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)
Videoman
22.02.2018 14:16То что вы описали очень похоже на вейвлет Хаара, нет? Каждый следующий шаг раскладывает сигнал на низкочастотную и высокочастотную составляющие.
Qvxb Автор
22.02.2018 14:19-1А что именно похоже? Посмотрел бегло про преобразование Хаара, оно вообще не про спектральный анализ. Теперь посмотрел подробнее, действительно есть упоминания про частотный фильтр, но там насколько я понял другая идея
Videoman
22.02.2018 14:24+1Хаар вообще-то математик и он много чего понапридумывал. При чем тут преобразование? Я ж вам дал ссылку на вейвлет — очень похоже на ваш способ.
Qvxb Автор
22.02.2018 14:38Посмотрел по ссылке пункт: преобразование для одномерного сигнала. Там исходный сигнал раскладывается на 2 других, но чтобы их получить нужно складывать/вычитать соседние точки. Речи про производные там нет.
T_Sun
22.02.2018 16:12Да, мне тоже вспомнилось дискретное вейвлет-преобразование Хаара. Спектр сигнала при этом вычисляется изменением ширины вейвлет-функции: 2, 4, 8 и т.д. Функция умножается на исходный сигнал, получается «послойное» разложение сигнала на частоты. Может быть не в чистом виде так, но очень похоже.
Qvxb Автор
22.02.2018 16:22По вашей ссылке есть код. Сравните с тем что в статье.
T_Sun
23.02.2018 07:55Код — всего лишь способ реализации идеи (алгоритма), одну идею можно реализовать разным кодом, нет смысла сравнивать коды напрямую.
Я же написал, что ваша идея мне напоминает вейвлет-преобразование и описал какое именно.
Я не утверждал, что у вас вейвлеты, всего лишь сказал, что ваша идея похожа на вейвлеты Хаара.
Так это или нет решать вам, вы лучше меня представляете вашу же идею выделения спектра и её реализацию в виде программы.
Sadler
22.02.2018 16:09+2Зачем вся эта цветомузыка? Вместо неё должен быть пример спектрограммы, полученной с помощью алгоритма. Можно в сравнении с классическим STFT.
Watcover3396
22.02.2018 17:37+2Мне показалось, или Вы пере-изобрели Преобразование Фурье?
GermanRLI
22.02.2018 19:10Автор по крайней мере попытался его переизобрести, но куда-то потерял фазу.
Действительно надо посмотреть на результат работы с одной/несколькими синусоидами.
Qvxb Автор
22.02.2018 17:40Здесь же Гиктаймс а не хабр, а эта цветомузыка и есть спектрограмма.(это ответ для Sadler)
Sadler
22.02.2018 19:30Претензия к Вам лишь в том, что вы выбрали крайне неинформативное представление результатов, и я не знаю, продиктовано это отсутствием опыта или желанием скрыть общую картину. На хабре Вы бы ещё добавили код этого безобразия.
LorDCA
23.02.2018 00:27Вообще, не думаю что человек воспринимает звук как разложение в спектр, скорее звук состоит из звуковых образов, которые состоят из волн разной длинны. Соответственно громкость какого либо звука это скорее средняя громкость всех волн из которых звук состоит. Но пока что не понятно какие параметры составляют звуковой образ, возможно средняя частота, стандартное отклонение или что то еще.
Человек воспринимает звук субъективно. То есть по сути мы понятия не имеем как конкретный человек слышит тоже само что и мы. Описывает он это теми же самыми словами, но что по факту, мы можем только гадать.
Громкость в ухе передается тактом от 20 до 100 гц.
А про звуковые образы и обобщение различных сенсоров (ака сенсорные нейроны, которых миллионы), можно почитать кучу книг, и написать еще не меньшую кучу. :)Qvxb Автор
23.02.2018 11:41Да, субьективно. Но все различают скачки (удары) и монотонный звук. И также различают гласные, согласные. Т.е мозг таки имеет(учится иметь) стандартные критерии в сигнале, на основании которых делает вывод что это за звук.
А про передачу громкости в ухе, имеется ввиду частота спайков нейронов?
aamonster
23.02.2018 07:26- Этап до получения промежуточный0… промежуточный8 очень долго описывается, а на самом деле это просто пачка ФНЧ. Если это понять — и проследить алгоритм (предсказав его поведение) проще, и ускорить можно.
- Как вам посоветовали, попробуйте на смеси двух синусоид. Сразу скажу, тут будет косяк: если обе частоты попадают в один слой — ваш детектор выдаст одну частоту (с переменной амплитудой).
- Хотелось бы видеть сравнение по быстродействию с БПФ.
В общем, для цветомузыки сойдёт (если не медленней бпф, тогда нет особого смысла), а вот для анализа сигналов я бы уже не рискнул применять.
Бонус (раз уж цветомузыка): фильтр типа y(t) = 0.9y(t-T)+0.1x(t) выделяет частоту 1/T и её гармоники, если сделать пачку таких фильтров для 12 нот какой-то октавы — можно отрисовать 12 красивых осциллограмм, дёргающихся в такт нотам.
Qvxb Автор
23.02.2018 11:361)Ну конечно это пачка фильтров. Но ведь все дело именно в том как фильтр работает.
2)Про смесь 2 синусоид, если сложить 2 синусоиды частоты которых отличаются меньше чем в 2 раза то никакой алгоритм не даст 2 этих синусоиды. Поправьте если не так.
По поводу фильтра типа y(t) = 0.9y(t-T)+0.1x(t), а есть примеры работы этого фильтра которые можно посмотреть?aamonster
23.02.2018 14:222 — преобразование Фурье даст (оно для этого и предназначено)
Пример — с ходу не найду, если покопаться в _очень старых завалах_ — возможно, накопаю исходник на паскале лохматого года под DOS, цифрующий звук с Sound Blaster Pro и отображающий 12 графиков в 13h видеорежиме VGA (320*240, 256 цветов :-))
Наверное, это уже даже не запустить без DosBox… Но визуально эффект красивый, нота в мелодии отображается сигналом на одном из графиков, причём форма сигнала — не синусоида (есть гармоники)
Да, регулировка «чувствительности» (добротности фильтров) — баланс между 0.9 и 0.1. Сейчас уже не вспомню — возможно, у меня было 0.95 и 0.05, или как-то так.Qvxb Автор
23.02.2018 15:562 — но только я сомневаюсь что оно даст 2 чистые синусоиды, скорее весь интервал между частотами будет не 0 (я пробовал юзать FFT раньше, так что примерно представляю какой оно выдает результат).
А как называется этот фильтр?aamonster
23.02.2018 16:19Чистые синусоиды (точнее, амплитуды/фазы двух синусоид) — даст. Единственное — если окно не кратно частоте синусоиды, она ещё на какие-то частоты залезает (см. «боковые лепестки», «растекание спектра»), лечится это с помощью оконной функции.
Фильтр — гребенчатый фильтр или что-то типа того. Крайне примитивная и быстро считающаяся вещь, но именно для красивого отображения нот хорош.
Dr_Dash
23.02.2018 16:55Однозначно, различит. Предел различимости частот 1/T, где Т — длительность сэмпла. В два раза выше должна быть частота взятия отсчётов, чем высшая измеряемая частота в спектре. А практически, для качественных измерений частота взятия отсчётов берётся в 3 раза выше. Я как практик говорю, сам разрабатываю.
Qvxb Автор
23.02.2018 17:36Да тут речь не про теорему Найквиста/Котельникова, а про результат после анализа. Т.е на мой взгляд после анализа суммы двух близких частот F1/F2<2, FFT даст весь набор между F1 и F2, возможно у F1 и F2 будут больше амплитуды, но будет и много лишнего.
Dr_Dash
23.02.2018 17:43нет, результат довольно чистый. Обязательно накладывается оконная функция. чтобы, как написали выше не было «растекания спектра». И результат настолько чистый, что побочные гармоники не наблюдаются практически.
Qvxb Автор
23.02.2018 18:45Если взять например, волну с одной единственной частотой, но периодически меняющейся амплитудой, то FFT дает вот что:
вот как выглядит сигнал:
Если конечно windows player корректно юзает FFTDr_Dash
24.02.2018 08:50Странно и нелепо. Вы бы написали простенькую программу да и анализированли. вот FFT отработало на массиве 16к точек, три синусоиды с частотами 1000, 1004, 1008
Даже окно не пришлось накладывать, но для реальных сигналов лучше наложить, использую окно Хэмминга если интересно.Qvxb Автор
24.02.2018 12:39Ведь этот алгоритм написан уже много раз, зачем писать еще 1? Можете дать ссылку на программу которая на ваш взгляд корректно отрабатывает FFT. Например вот эта программа норм: ua3vvm.qrz.ru/gram.htm?
Dr_Dash
24.02.2018 12:52Можете дать ссылку на программу которая на ваш взгляд корректно отрабатывает 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); }Qvxb Автор
24.02.2018 13:14Ок, спасибо, на самом деле результат на Вашем скриншоте выглядит круто. Повожусь с кодом который Вы дали.
Qvxb Автор
24.02.2018 13:50int 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 что в массиве IDr_Dash
24.02.2018 15:28нужно получать DP__Amplitude_FFT, и там уже спектр будет
Qvxb Автор
24.02.2018 15:47пробовал и с 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]);
Qvxb Автор
24.02.2018 15:57может я неправильно обозначил u4x, я так понял что это целый тип, обозначил его как short. #define u4x short
Dr_Dash
24.02.2018 18:23Большие числа у вас потому что большие амплитуды синусов, оставьте 1.
Qvxb Автор
24.02.2018 19:39Поставил амплитуду 1, в результате получилось куча значений больше 1, даже есть пара значений больше 1000. Правильно ли я понимаю что в массиве Desination[i], амплитуда гармоники — значение Desination[i], а частота — индекс i?
Dr_Dash
24.02.2018 21:04Попробуйте одну синусоиду,
Signal[i]= sin(2.0*3.1415*float(i)*1000.0/16384);
и почему вы делаете так
if(Desination[i]>1)?
про Desination[i] понимаете верноQvxb Автор
25.02.2018 12:05Все, заработало. Потестил разные варианты, работает очень даже норм.
Desination[i]>1, это чтобы не выводить 0 значения.
aamonster
23.02.2018 16:21s/240/200/
(на случай, если придёт кто-то, кто помнит историю, и заметит мою ошибку)
Borjomy
25.02.2018 11:23И что касается фильтров и что касается преобразований Фурье, все они анализируют определенный временной участок. Чем более крутой фильтр, выше его порядок, тем дольше он «устанавливается». То-же относится и к Фурье. Преобразование Фурье дает среднеее значение амплитуды и фазы анализируемого фрагмента. Чем больший участок времени вы берете, тем точнее разделяются частоты, но за все надо платить, и вычисленные вариации амплитуды сглаживаются. Вообще, преобразование Фурье применяется для анализа бесконечной волновой функции, речь и музыка к таковым не относятся. Выше кто-то исследовал работу Фурье при изменении амплитуды. Господа, вы забываете, что изменение амплитуды, т.е амплитудная модуляция, приводит к изменению амплитуд и соседних частот.
Если вы берете какой-то фрагмент для анализа, то в него попадает один период самой низкой частоты и 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
И не забывайте накладывать окна (Хемминга например). Преобразования Фурье очень чувствительны к выбросам на концах анализируемого массива.
Dr_Dash
Вы не пробовали 3-4 заранее известных синусоиды смиксовать и подать и посмотреть что получится? Очень интересно, если можете, добавьте в статью.
Qvxb Автор
ок, добавлю