Хабр Курсы для всех
РЕКЛАМА
Практикум, Хекслет, SkyPro, авторские курсы — собрали всех и попросили скидки. Осталось выбрать!
Хабр доступен 24/7 благодаря поддержке друзей
Комментарии 42
Вы не пробовали 3-4 заранее известных синусоиды смиксовать и подать и посмотреть что получится? Очень интересно, если можете, добавьте в статью.
То что вы описали очень похоже на вейвлет Хаара, нет? Каждый следующий шаг раскладывает сигнал на низкочастотную и высокочастотную составляющие.
А что именно похоже? Посмотрел бегло про преобразование Хаара, оно вообще не про спектральный анализ. Теперь посмотрел подробнее, действительно есть упоминания про частотный фильтр, но там насколько я понял другая идея
Хаар вообще-то математик и он много чего понапридумывал. При чем тут преобразование? Я ж вам дал ссылку на вейвлет — очень похоже на ваш способ.
Посмотрел по ссылке пункт: преобразование для одномерного сигнала. Там исходный сигнал раскладывается на 2 других, но чтобы их получить нужно складывать/вычитать соседние точки. Речи про производные там нет.
Да, мне тоже вспомнилось дискретное вейвлет-преобразование Хаара. Спектр сигнала при этом вычисляется изменением ширины вейвлет-функции: 2, 4, 8 и т.д. Функция умножается на исходный сигнал, получается «послойное» разложение сигнала на частоты. Может быть не в чистом виде так, но очень похоже.
По вашей ссылке есть код. Сравните с тем что в статье.
Код — всего лишь способ реализации идеи (алгоритма), одну идею можно реализовать разным кодом, нет смысла сравнивать коды напрямую.
Я же написал, что ваша идея мне напоминает вейвлет-преобразование и описал какое именно.
Я не утверждал, что у вас вейвлеты, всего лишь сказал, что ваша идея похожа на вейвлеты Хаара.
Так это или нет решать вам, вы лучше меня представляете вашу же идею выделения спектра и её реализацию в виде программы.
Я же написал, что ваша идея мне напоминает вейвлет-преобразование и описал какое именно.
Я не утверждал, что у вас вейвлеты, всего лишь сказал, что ваша идея похожа на вейвлеты Хаара.
Так это или нет решать вам, вы лучше меня представляете вашу же идею выделения спектра и её реализацию в виде программы.
Зачем вся эта цветомузыка? Вместо неё должен быть пример спектрограммы, полученной с помощью алгоритма. Можно в сравнении с классическим STFT.
Мне показалось, или Вы пере-изобрели Преобразование Фурье?
Здесь же Гиктаймс а не хабр, а эта цветомузыка и есть спектрограмма.(это ответ для Sadler)
Вообще, не думаю что человек воспринимает звук как разложение в спектр, скорее звук состоит из звуковых образов, которые состоят из волн разной длинны. Соответственно громкость какого либо звука это скорее средняя громкость всех волн из которых звук состоит. Но пока что не понятно какие параметры составляют звуковой образ, возможно средняя частота, стандартное отклонение или что то еще.
Человек воспринимает звук субъективно. То есть по сути мы понятия не имеем как конкретный человек слышит тоже само что и мы. Описывает он это теми же самыми словами, но что по факту, мы можем только гадать.
Громкость в ухе передается тактом от 20 до 100 гц.
А про звуковые образы и обобщение различных сенсоров (ака сенсорные нейроны, которых миллионы), можно почитать кучу книг, и написать еще не меньшую кучу. :)
- Этап до получения промежуточный0… промежуточный8 очень долго описывается, а на самом деле это просто пачка ФНЧ. Если это понять — и проследить алгоритм (предсказав его поведение) проще, и ускорить можно.
- Как вам посоветовали, попробуйте на смеси двух синусоид. Сразу скажу, тут будет косяк: если обе частоты попадают в один слой — ваш детектор выдаст одну частоту (с переменной амплитудой).
- Хотелось бы видеть сравнение по быстродействию с БПФ.
В общем, для цветомузыки сойдёт (если не медленней бпф, тогда нет особого смысла), а вот для анализа сигналов я бы уже не рискнул применять.
Бонус (раз уж цветомузыка): фильтр типа y(t) = 0.9y(t-T)+0.1x(t) выделяет частоту 1/T и её гармоники, если сделать пачку таких фильтров для 12 нот какой-то октавы — можно отрисовать 12 красивых осциллограмм, дёргающихся в такт нотам.
1)Ну конечно это пачка фильтров. Но ведь все дело именно в том как фильтр работает.
2)Про смесь 2 синусоид, если сложить 2 синусоиды частоты которых отличаются меньше чем в 2 раза то никакой алгоритм не даст 2 этих синусоиды. Поправьте если не так.
По поводу фильтра типа y(t) = 0.9y(t-T)+0.1x(t), а есть примеры работы этого фильтра которые можно посмотреть?
2)Про смесь 2 синусоид, если сложить 2 синусоиды частоты которых отличаются меньше чем в 2 раза то никакой алгоритм не даст 2 этих синусоиды. Поправьте если не так.
По поводу фильтра типа y(t) = 0.9y(t-T)+0.1x(t), а есть примеры работы этого фильтра которые можно посмотреть?
2 — преобразование Фурье даст (оно для этого и предназначено)
Пример — с ходу не найду, если покопаться в _очень старых завалах_ — возможно, накопаю исходник на паскале лохматого года под DOS, цифрующий звук с Sound Blaster Pro и отображающий 12 графиков в 13h видеорежиме VGA (320*240, 256 цветов :-))
Наверное, это уже даже не запустить без DosBox… Но визуально эффект красивый, нота в мелодии отображается сигналом на одном из графиков, причём форма сигнала — не синусоида (есть гармоники)
Да, регулировка «чувствительности» (добротности фильтров) — баланс между 0.9 и 0.1. Сейчас уже не вспомню — возможно, у меня было 0.95 и 0.05, или как-то так.
Пример — с ходу не найду, если покопаться в _очень старых завалах_ — возможно, накопаю исходник на паскале лохматого года под DOS, цифрующий звук с Sound Blaster Pro и отображающий 12 графиков в 13h видеорежиме VGA (320*240, 256 цветов :-))
Наверное, это уже даже не запустить без DosBox… Но визуально эффект красивый, нота в мелодии отображается сигналом на одном из графиков, причём форма сигнала — не синусоида (есть гармоники)
Да, регулировка «чувствительности» (добротности фильтров) — баланс между 0.9 и 0.1. Сейчас уже не вспомню — возможно, у меня было 0.95 и 0.05, или как-то так.
2 — но только я сомневаюсь что оно даст 2 чистые синусоиды, скорее весь интервал между частотами будет не 0 (я пробовал юзать FFT раньше, так что примерно представляю какой оно выдает результат).
А как называется этот фильтр?
А как называется этот фильтр?
Чистые синусоиды (точнее, амплитуды/фазы двух синусоид) — даст. Единственное — если окно не кратно частоте синусоиды, она ещё на какие-то частоты залезает (см. «боковые лепестки», «растекание спектра»), лечится это с помощью оконной функции.
Фильтр — гребенчатый фильтр или что-то типа того. Крайне примитивная и быстро считающаяся вещь, но именно для красивого отображения нот хорош.
Фильтр — гребенчатый фильтр или что-то типа того. Крайне примитивная и быстро считающаяся вещь, но именно для красивого отображения нот хорош.
Однозначно, различит. Предел различимости частот 1/T, где Т — длительность сэмпла. В два раза выше должна быть частота взятия отсчётов, чем высшая измеряемая частота в спектре. А практически, для качественных измерений частота взятия отсчётов берётся в 3 раза выше. Я как практик говорю, сам разрабатываю.
Да тут речь не про теорему Найквиста/Котельникова, а про результат после анализа. Т.е на мой взгляд после анализа суммы двух близких частот F1/F2<2, FFT даст весь набор между F1 и F2, возможно у F1 и F2 будут больше амплитуды, но будет и много лишнего.
нет, результат довольно чистый. Обязательно накладывается оконная функция. чтобы, как написали выше не было «растекания спектра». И результат настолько чистый, что побочные гармоники не наблюдаются практически.
Если взять например, волну с одной единственной частотой, но периодически меняющейся амплитудой, то FFT дает вот что:
вот как выглядит сигнал:
Если конечно windows player корректно юзает FFT
вот как выглядит сигнал:
Если конечно windows player корректно юзает FFT
Странно и нелепо. Вы бы написали простенькую программу да и анализированли. вот FFT отработало на массиве 16к точек, три синусоиды с частотами 1000, 1004, 1008
Даже окно не пришлось накладывать, но для реальных сигналов лучше наложить, использую окно Хэмминга если интересно.
Даже окно не пришлось накладывать, но для реальных сигналов лучше наложить, использую окно Хэмминга если интересно.
Ведь этот алгоритм написан уже много раз, зачем писать еще 1? Можете дать ссылку на программу которая на ваш взгляд корректно отрабатывает FFT. Например вот эта программа норм: ua3vvm.qrz.ru/gram.htm?
Можете дать ссылку на программу которая на ваш взгляд корректно отрабатывает 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 будут больше амплитуды, но будет и много лишнего.
Ок, спасибо, на самом деле результат на Вашем скриншоте выглядит круто. Повожусь с кодом который Вы дали.
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
нужно получать DP__Amplitude_FFT, и там уже спектр будет
пробовал и с 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]);может я неправильно обозначил u4x, я так понял что это целый тип, обозначил его как short. #define u4x short
s/240/200/
(на случай, если придёт кто-то, кто помнит историю, и заметит мою ошибку)
(на случай, если придёт кто-то, кто помнит историю, и заметит мою ошибку)
Скажите, а вы не пробовали замерять скорость работы вашего алгоритма и быстрых преобразований Фурье? А также сравнить результаты работы вашего алгоритма и БПФ?
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Еще один способ разложения сигнала в спектр