Поиск:

Ответ в темуСоздание новой темы Создание опроса
> Как определить точную частоту звукового сигнала, какой алгоритм лучше использовать ? 
:(
    Опции темы
Santik
Дата 20.3.2012, 06:47 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



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

Работаем с модулем FFT !!!
Можно 2 раза - т.к. одну точку ты уже нашёл через FFT.  Если алгоритм расчёта параболы рассчитан на равномерную сетку, тогда три раза интеграл считать.
Но это лучше перебора с шагом 0.5.  И лучше с 44100 работать.
На увеличение быстродействия особо не надейся  smile 
Прикинь сколько вообще ты можешь итераций себе позволить.
Учти, что на коротком сигнале ты не получишь хорошую точность! И ничего не поможет! 
Я думаю точность 0.1 Гц для 6 струны получишь при длительности сигнала 1 секунда.
PM MAIL   Вверх
Pavia
Дата 20.3.2012, 19:02 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 418
Регистрация: 6.12.2008

Репутация: 11
Всего: 12



Самый простой способ увеличить в 10-20 раз это взять библиотеку FFTW, она быстрее обычной самописной реализации.
Берем 2 буфера
1)  N=2048 отсчётов
2)  NewN=64*1024 отсчётов
В первый пишем сигнал пишем сигнал потом умножаем на оконную функцию лучше Ланцоша (sinc)
Копируем результат во второй остававшийся часть заполняем нулями.  
После FFT и ищем максимум.
i*Fd/(NewN); частоту определяем как.
Точность будет  0.1 Гц.
Больше ничего не требуется. Без FFTW по скорости не отличается от предыдущего алгоритма.


Цитата
 переделки алгоритма 
 Извени, но у меня есть и другие вещи, которыми я могу заниматься. 
PM MAIL   Вверх
Santik
Дата 21.3.2012, 14:39 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Сейчас проверил ещё один интересный способ уточнения частоты.
Алгоритм такой:
1. FFT и находим основную частоту
2. Определяем сигнал с линейно меняющейся частотой  Fmax +- 0.5 Гц    (44100)
3. FFT  этого сигнала
4. Перемножаем вектора
5. Реверсное FFT--> X(i)
5. Отношение Х(N/2)/X(1) даст поправку к определённой в п.1 частоте

Проверил на тестовом сигнале. Точность 0.001 Гц !!!  (Fd=44100)

PM MAIL   Вверх
Kpeved
Дата 21.3.2012, 22:48 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 42
Регистрация: 26.2.2012

Репутация: нет
Всего: нет



Santik, а как по производительности? Сколько итераций это все приблизительно займет ?  И что значит "Сигнал  с линейно меняющейся частотой"? 
PM MAIL   Вверх
Santik
Дата 21.3.2012, 23:40 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



А итераций нет совсем. Грубо говоря - всё забито в прямое и обратное преобразование Фурье.
Алгоритм я уже приводил на Фортране. Там нужно только ввести сигнал с переменной частотой: f=f0+(0.5)/T*t
Все примеры для 44100, Т=1сек.
Пусть мы из FFT определили  fmax= 328 Гц . Если значение при 329 больше чем в 327 то нам надо двигаться в +.
f(t)=328+(0.5)/T*t , т.е. частота меняется линейно во времени от 328 до 328.5 Гц
Интеграл от 2*Pi*f(t) даст фазу сигнала Ф(t). 
 cos(Ф(t)) +i sin(Ф(t)) - это и будет наш сигнал с линейно меняющейся частотой.  
Далее делаем свёртку этого сигнала с исходным сигналом (см. программу)
(умножаем вектора результатов FFT  и делаем обратное FFT
Максимум модуля свёртки (или точная частота) гарантированно будет в этом диапазоне.

Завтра приведу данные по производительности. Программу полноценную мне конечно лень писать. Но для примера текст отладочной программы приведу, чтобы можно разобраться было.
Ну а дальше - сам. smile 
PM MAIL   Вверх
Santik
Дата 22.3.2012, 08:31 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Данные по быстродействию:
Fd=44100
N=44100  (1 сек)
Тестовый сигнал  329.123456 Гц с шумом.
600 значений частот рассчитаны за 53 секунды. (P-4 3ГГц) Погрешность 0.001 Гц.
PM MAIL   Вверх
Santik
Дата 22.3.2012, 18:49 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Код

integer*4 function GSM19_3(Wav_name)
use msfwina
use IMSL
USE DFWIN
USE DFPORT 

Type(T_COORD) wpos
character*100    Wav_name
integer*4 i,imax,j,k,N,N2,Ng
double complex , dimension (440000) :: X,FX, X1, X2
real*8  N1,xr,xi,f0,fk,fmax,df,f10,f0x,Fd,Fd1,wt,wt0,eR,ey,t
real*8  h,x0,y_1,y0,y1
integer*2 , dimension (400000) :: Vib1
Real*8 , dimension (40) :: A
integer*2 xxx
double complex  ez

OPEN (8, FILE = '44k-1.txt'C, ACTION='READ',FORM='UNFORMATTED', ACCESS = 'DIRECT',RECL=2)
OPEN (7, FILE = '1004.txt'C, ACCESS = 'SEQUENTIAL' , STATUS = 'REPLACE')

Fd=44100  !частота дискретизации
Fd1=1/Fd
N=44100*1!                44100 - 1 секунда
N1=1./N
N2=N/2
do  i=1, N !  Cчитываем данные 
read(8,REC=i+44100+12900+30000 )xxx !+44100+12900+30000
Vib1(i)=(xxx)
er=real(Vib1(i))
X(i)= cmplx(er,0.0)
end do

ret = MessageBox (hwndDlg,'START!', "ВНИМАНИЕ!"C, MB_OK) 


go to 200            ! Обход тест-сигнала
f0x=298.12345678  !     +(k-1)/1200.  
ar= 1+0.5*rand (0)*(-1)**k
fi=2*3.1415926*rand (0)*(-1)**k
do  i=1, N 
t=(i-1)*Fd1
xr=2*3.1415926*f0x*t+ fi
er= ar*(cos(xr)+sin(xr)) + 0.01*rand (0)*(-1)**k
X(i)= cmplx(er,0.0)
end do
200 continue

    call DFFTCF(N,X,FX)  ! БПФ
! Здесь должна быть подпрограммы поиска мах модуля FX
!Пока f0= 328
! Максимум в 328 Гц и уточняем "Вверх"

f0= 328.0      !начальная частота
df=0.5

do i=1,N    !Формируем сигнал  свипа 
t=(i-1)*Fd1
f10= 2*3.1415926*(f0 + df/(2*N)*Fd*t)*t   !Фаза сигнала
    xr=  cos(f10) 
    xi=  sin(f10)
X(i)= cmplx(xr,xi)
end do

       call DFFTCF(N,X,X1)        !БПФ свипа
  
 do i=1,N
 X2(i)= X1(i)*DCONJG (FX(i))    !  DCONJG  - комплексное сопряжение
end do 

    call DFFTCB(N,X2,X1)  !  в Re(X1) - свёртка
    
    er=CDABS(X1(N/2))/CDABS(X1(1))
df=-0.861*er+0.85504     !Находим поправку

f11=f0+df !  Истинная частота

        write(7,1032)k,df,df1,f0x, f11, f0x-f11,f0x-f0,er
1032        format(i6,6F15.7,E12.5)

    close(7)
    close(8)
gsm19=N

return
end

PM MAIL   Вверх
Santik
Дата 25.3.2012, 06:31 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Вот ещё ссылка по точности:
http://dsplib.ru/forum/viewtopic.php?f=7&a...42&start=10
PM MAIL   Вверх
Kpeved
Дата 26.3.2012, 00:36 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 42
Регистрация: 26.2.2012

Репутация: нет
Всего: нет



Наконецто набрася храбрости прочитать код на ещё одном неизвестном мне языке .) По коду мне все ещё непонятно следующее :
Цитата

67:  X2(i)= X1(i)*DCONJG (FX(i))    !  DCONJG  - комплексное сопряжение 
 - что такое FX(i) и вообще что оно делает ?   X1(i)- как я понял структура с реальной и мнимой частью ? Или модуль и фаза ?
Ещё не знаю - есть ли в фортране перегрузка операторов - но как именно он умножает операнды ? 
Да, и что делает DCONJG тоже не ясно )
Цитата

72:    er=CDABS(X1(N/2))/CDABS(X1(1))
73: df=-0.861*er+0.85504     !Находим поправку
 Что делает CDABS  ? И что за магические числа ?)
Вроде бы по коду все .

Ещё по быстродействию - Если мы будем использовать ЧД 8К при 1024 значениях - будет ли у нас нужная точность ? Получается что на 1 отсчет придется по 8Гц и +- надо брать намного больше .
PM MAIL   Вверх
Santik
Дата 26.3.2012, 06:53 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Цитата(Kpeved @ 26.3.2012,  00:36)
Наконецто набрася храбрости прочитать код на ещё одном неизвестном мне языке .) 

Даже это уже положительный момент!!!
Да и что может быть проще Фортрана?
Приивлекает богатейшая библиотека IMSL.
 БПФ, к примеру, - с любым количеством точек и очень быстро!
 DCONJG  -функция комплексного сопряжения (комплексного вектора FX -  в данном случае в FX  - БПФ исходного сигнала ) 
CDABS - -функция  возвращае модуль  комплексного числа (D- двойная точность)

По быстродействию сказать ничего не могу.
При изменении N "магические числа" могут быть другими!
По точности см. 
http://dsplib.ru/forum/viewtopic.php?f=7&a...42&start=10 
Выбираешь нужную погрешность и находишь М и потом N - число точек (для 6 струны - обязательно, т.к. НЧ хуже по точности)

Я пока этот алгоритм не "обкатыва". Надеюсь совместными усилиями ....

Это сообщение отредактировал(а) Santik - 26.3.2012, 07:39
PM MAIL   Вверх
Kpeved
Дата 26.3.2012, 13:31 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 42
Регистрация: 26.2.2012

Репутация: нет
Всего: нет



Santik, ну а что вообще за числа в 73 ? Как их рассчитывать ?
Как я понял М и N рассчитываются по последней формуле ?  Тогда если Т=0.128 ,f= 70Гц, то M =8.96. Погрешность около 0.1 Гц .  И ??)
PM MAIL   Вверх
Santik
Дата 26.3.2012, 15:05 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Цитата(Kpeved @ 26.3.2012,  13:31)
Santik, ну а что вообще за числа в 73 ? Как их рассчитывать ?
Как я понял М и N рассчитываются по последней формуле ?  Тогда если Т=0.128 ,f= 70Гц, то M =8.96. Погрешность около 0.1 Гц .  И ??)

Определение коэффициентов: 
Снимаем характеристику  er=CDABS(X1(N/2))/CDABS(X1(1)) при расстройке частоты от 0 до 0.5 с шагом 0.01 Гц.
Находим ошибку. Находим тренд функции Ошибка(er).  (я для начала в Excel  делал)  Получаем df=-0.861*er+0.85504     
По полученному значению er находим ошибку df.
У меня сейчас:
df=0.11356244E+01*er**7-0.56608063E+01*er**6+0.11240384E+02*er**5-0.11300708E+02*er**4+ 0.64610102E+01*er**3-0.25383018E+01*er**2  -0.60435858E-01*er+ 0.72327669E+00  smile 

...  И ??)

Это говорит о том, что при данной частоте и времени сигнала ты определишь частоту с погрешностью не менее 0.1 Гц .
Если тебе нужна меньшая погрешность - придется увеличивать Т  и другого выхода нет!
PM MAIL   Вверх
Santik
Дата 26.3.2012, 16:40 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Для Fd=8000 Гц  и  N=1024

df=-0.69019186E+01*er+ 0.68520498E+01

При частоте 78 Гц метод даёт погрешность 0.2 Гц
PM MAIL   Вверх
Santik
Дата 26.3.2012, 20:16 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



В файле некоторые результаты вычисления при Fd=8000 Гц.

Присоединённый файл ( Кол-во скачиваний: 10 )
Присоединённый файл  8000_Гц.rar 33,87 Kb
PM MAIL   Вверх
Santik
Дата 29.3.2012, 11:12 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 45
Регистрация: 13.3.2012
Где: Мирный (Якутия)

Репутация: нет
Всего: нет



Исследуем полифонический сигнал (файл 44k-all.txt)!

Это сообщение отредактировал(а) Santik - 29.3.2012, 19:09

Присоединённый файл ( Кол-во скачиваний: 10 )
Присоединённый файл  1029_3.rar 9,94 Kb
PM MAIL   Вверх
Страницы: (9) Все « Первая ... 5 6 [7] 8 9 
Ответ в темуСоздание новой темы Создание опроса
Правила форума "Алгоритмы"

maxim1000

Форум "Алгоритмы" предназначен для обсуждения вопросов, связанных только с алгоритмами и структурами данных, без привязки к конкретному языку программирования и/или программному продукту.


Если Вам понравилась атмосфера форума, заходите к нам чаще! С уважением, maxim1000.

 
1 Пользователей читают эту тему (1 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема | Алгоритмы | Следующая тема »


 




[ Время генерации скрипта: 0.0935 ]   [ Использовано запросов: 20 ]   [ GZIP включён ]


Реклама на сайте     Информационное спонсорство

 
По вопросам размещения рекламы пишите на vladimir(sobaka)vingrad.ru
Отказ от ответственности     Powered by Invision Power Board(R) 1.3 © 2003  IPS, Inc.