Мне кажется FFT можно для синусов и косинусов можно по простой формуле посчитать.
Вот я для 7 гармоник прикинул (в файле).
Я не понял, причем сдесь Ng=1???
FFT вне этого цикла !
Если FFT  не делать - в цикле останется умножение 2-х векторов и реверсное FFT .



Добавлено @ 10:48
Файл

Это сообщение отредактировал(а) Santik - 17.3.2012, 10:52

Присоединённый файл ( Кол-во скачиваний: 2 )
 1007_1.xlsx 838,35 Kb

Дыра в быстродействии не в ФФТ а в синусах и косинусах  . При нахождении 1го пика мы мы M-раз , в зависимости от точности , проходим по всему массиву длинной N и каждый раз считаем син и кос . Потом ещё берём корень...
По быстродействию должен получится как ФФТ , но получается намного медленнее .

Может можно применить формулу эйлера или ещё что то ?

Неее... Ты не прав!  Просто каждый раз приходится вычислять "длинный" интеграл внутри которого умножения на син/кос.
А то что М-проходов - сделай итерационную процедуру поиска экстремума. У меня до 0.001 за 5-6 итераций сходилось...  А так "тупо" никто не делает. Pavia  это же для примера написал...

Добавлено через 13 минут и 24 секунды
И мне кажется надо попробывать как vedun предлагал, но только использовать все гармоники. Только полученные не забывай поделить на номер. И среднюю взять. Наверное точности хватит.

 А как вы собираетесь использовать все гармоники?

Ну первую мы вычислили... А дальше "тупо" умножаем на 2 , 3, .... 11 . И методом параболы уточняем.
Чем выше гармоника, тем точнее расчёт...

Это сообщение отредактировал(а) Santik - 17.3.2012, 12:02

 Но ведь эти гармоники в общем случае не кратны друг другу. Как вы собираетесь это учесть ?

Santik, 
Вы считаете через FFT свёртку. Для sin и cos свёртку линейная.

Не кажется так и есть. Вот корреляцию для вашего случая можно считать со сложностью N*Ng.
Но это не даёт преимуществ. Надо выносить FFT из цикла по k.


Santik, 
Я интерационно не считал, так как с полосой был не уверен. Сейчас уверен, да можно интерационо. Только вот полоса пропускания 1 Гц поэтому интонационно можно только когда у нас такая точность 1 Гц.  В виду того что у нас массив менее секунды точность не дотягивает до 1 Гц.
Так что тут прирост будет в 2 раза. 

А вот свёртку можно ускорить в виду того что у меня используется FT где sin и cos вычисляется на каждом шагу.
Есть алгоритм Герцеля. В виду того что он плохо описан в литературе, дальнейшее что я опишу  может отличаться или не соответствовать.
Так вот для быстрого расчёта можно использовать свойство косинус суммы и синус суммы. 



Здесь Cos(-2*Pi*i*w*N1); можно представить в виде cos(-2*Pi*(i-1)*N1+-2*Pi*1*N1)=переобозначим=cos(a_I_1+b)=cos(a_I_1)*cos(b)-sin(a_I_1)*sin(b)
sin(a_I_1+b)=sin(a_I_1)*cos(b)+cos(a_I_1)*sin(b)

1)sin(b)
cos(b) 
Вычисляются вне цикла.
2)
sin(a_I_1)
cos(a_I_1)
Известно с предыдущей итерации цикла. Раз в 10 должно ускориться.

Добавлено через 12 минут и 34 секунды
А вообще всё это ерунда, берём FFT и немного переделываем поделив  начальную частоту в N раз.
Тем самым мы повысим точность и скорость.

А почему не кратны? Это же колебания струны. Обязаны быть кратными!
Я проверял. Правда если вы не на 44100 Гц  смотрели , то возможно, что пролезают гармоники с частотами больше Найквиста, которые благополучно в низкочастотную отражаются...

Это сообщение отредактировал(а) Santik - 17.3.2012, 15:55

Так вот я  и хочу одно FFT из цикла убрать. Останется в цикле перемножение векторов и обратное FFT.
Тоже конечно не очень...
Так где мне посмотреть это "так и есть"??? У меня вылезает что-то типа SIN(X)/X.
Велосипед изобретать не хочется.

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

Исходя из задачи  .  Не думаю, что при таком уровне измерений мы дисперсию увидим.  Надо считать.

vedun,  докладываю: Дисперсии не обнаружено!



А если и есть, то совсем маленькая и почти совсем незаметная...  

Это сообщение отредактировал(а) Santik - 17.3.2012, 17:53

Я что-то пропустил, наверное! Новые методы появились???  

Pavia, а можно поподробнее как именно нужно переделать ? И зачем нам частоту делить ?

Pavia похоже погорячился...
Я вот испытал модифицированный метод Vedunа.
Через 3 точки Х(f0),Х(f1),Х(f2) (из FFT,  при f1- макс. значение) проводим параболу. Находим f3.
Далее модификация - считаем наш интеграл и находим Х(f3)
Ещё раз считаем интеграл Х(f3 +- h)      (+ или минус в зависимости от значений Х(f1) и Х(f3).
Через 3 точки снова проводим параболу. Находим искомую частоту.
Ошибка менене 0.1 Гц  на частоте 380 Гц. (Fd=44100)
Таким образом "длинный" интеграл считаем всего 2 раза!
Можно, конечно, еще раз эту процедуру повторить, чтобы получить ещё меньшую погрешность.


Это сообщение отредактировал(а) Santik - 19.3.2012, 21:56

Santik, а на чем проверяли ? На генераторе или реальном звуке ?  Мне кажется что там вообще не парабола - фронт и спады разные и она там никак не получится . 
И у вас как я понял 3 раза интеграл считается - для X(f3) , X(f3-h) , X(f3+h)  - потом находится методом параболы ? Т.есть мы просто считаем значения амплитуды в этих точках и строим параболу ? Не совсем точный метод  . Тем более если учесть искажённые сигналы с перегрузом или дисторшном или с чем нибудь ещё то там будет точно не парабола . 
Я все ещё надеюсь что  Pavia что нибудь напишет по поводу переделки алгоритма )