Programmist, сподобился я как бы ник свой оправдать - составил и испытал в маткаде алгоритм взятия табличных значений синуса :) Вроде работает. Для простоты в таблицу включил полный период. Значения выбираются с линейной интерполяцией. При длине таблицы M=1024 ошибка вышла не хуже чем в пятом знаке после запятой. Длина выходного массива N может быть любой. Если Вас это заинтересует, то txt-таблицы синуса (и косинуса) любой нужной длины могу выложить. Если в чём-то дал маху, будем надеяться, опытные коллеги вразумят :)

Пишу индексы вверху в скобках, т.к. не знаю, как тут на форуме нижние индексы впечатать.

Дано:
Массив чисел: отсчёты сигнала S⁽ⁿ⁾, где n = 0,1, ... N-1.
Количество отсчётов сигнала N и частота дискретизации Fdiscr.
Массив табличных значений синуса: SinTab^(m) = sin(m*(2пи / M)), где m = 0,1, ... M.

Найти:
Массив произведений S⁽ⁿ⁾*sin(2пи*F*t⁽ⁿ⁾), где
F - заданная "частота гетеродина" в Гц, t⁽ⁿ⁾ = n / Fdiscr есть время в секундах.

Введём в рассмотрение "нормированную частоту гетеродина" G = F / Fdiscr. Тогда фаза (весь аргумент) синуса запишется в виде:
2пи*F*t⁽ⁿ⁾ = n*2пи*G =n*(M*G)*(2пи / M).

Обозначим D = M*G. Видно, что n*D есть "эффективный номер m" для взятия значений синуса из нашей таблицы, а D есть его приращение. Чтобы из эфф. номера получать целое m в допустимых пределах 0,1, ... M, надо по мере роста n отбавлять из n*D целое число раз лишние М и отбрасывать дробную часть - такова главная идея.

Будем это делать в цикле, а чтобы уменьшить количество умножений, введём аккумулятор А для накапливания n*D сложением, и введём аккумулятор B для подсчёта лишних M. Выгодно задавать M в виде степеней двойки, тогда доступно быстрое умножение на М. Итак, вот весь алггоритм; после косой двойной скобки - комментарии:

// Определяем значения констант
G := F / Fdiscr
D := M*G
A := 0
B := 0

Для n от 0 до N-1 с шагом 1
// Тело цикла
C := A - M*trunc(B) // Тут у нас trunc(B) означает взятие целой части числа B
m := trunc(C)
A := A + D
B := B + G
S1 := SinTab^(m) // Берём из таблицы два соседних значения
S2 := SinTab^(m+1) // для линейной интерполяции
OutSin⁽ⁿ⁾ := S1 + (S2 - S1)*(C - m) // Это и есть наш sin(2пи*F*t⁽ⁿ⁾)
OutS⁽ⁿ⁾ := S⁽ⁿ⁾*OutSin⁽ⁿ⁾ // Это выходной сигнал
// Конец тела цикла

Извиняюсь за оффтоп, хотя... синус - он как бы тоже может быть типо "особенностью звучания", в общем достойный объект изучения... :)

Sinus Если в чём-то дал маху Имхо, особых "Махов" :-) это ж скорость звука, нет. :) Наоборот Programmist вооружён до зубов, что называется перо в руки.

з.ы. Или я попутал с теоремой чи постулатом Маха, чёт знакомое этот Мах, ну ладно, эт так к слову.

зызы К стати, к стати есть ошибка в алгоритме, попробуйте сами найти, если не выйдет, то подскажу. Ошибка типовая, её сразу видно. Но думки думками, а дела делами, вдруг вы её обошли или где предусмотрели неявно подумал я. Вот, ну и так как не нашёл таких неявных обходов, то сосбсно и прочекил, конешно никуда эта ошибка не делась, а сидит и мерзко ждёт своего часа, а он наступит рано или поздно. :-)

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



Учу ее читать и конвертировать файлы любого формата в wav и обратно. На спектре текстовый файл :) Может быть оно и не особо надо, но решил оставить эту возможность. Предстоит еще много работы с линейками и масштабированием, дальше пойдет полная оптимизация кода.

з.ы.
Имхо, здесь должно что-то получиться. Работа со статикой, сильно отличается от работы в "реальном" времени, с непонятным устройством, которое выдает неизвестно что. Любую ошибку можно поймать и исправить.
Те два спектра, которые выше, звуковая карта съест с аппетитом, да еще и своего добавит, а там потом пойди разберись, что было в начале, курица или яйцо :)

Mesh

Ошибка типовая, её сразу видно.

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

Вот думаю ещё подводные камни могут быть в использовании для аккумулятора плавающей точки хотя х. з. ИМХО лучше взять один аккумулятор целочисленный 32 бит или 64 бит и адресацию таблицы и дробную часть для управления интерполятором с него получать.

petr0v
С целыми числами конечно проще, тем более что они есть большие, 64 бита, но возникают непредвиденные ситуации, с преобразованием массивов, а время на эту операцию уходит больше, чем работа на прямую с Float.
Надо делать циклы и смотреть, что быстрей, без проверки ничего не понятно. Если идти совсем верной дорогой, то только на asm-e, принудительно включив SSE, 3D, или что там еще в этих камнях имеется, сейчас это не реально.
Вообще шикарно будет под это дело графические ядра задействовать, как Mesh советовал, да далеко мне пока до них, как до луны.

Programmist

Где-то вы выше писали что с float медленнее получается. При работе с float нужно учитывать потерю точности и т. п., ведь это напрямую может повлиять на точность генератора.

Не знаю архитектру процов современных, что будет происходить если например к целому положительному 64 бит числу 2^64-2 прибавить 5, какой результат получится?

Ошибка типовая, её сразу видно.

Спасибо. Но не уверен, что всё вижу. Вижу только, что ячейки A и B рано или поздно могут переполниться. Просто не хотел применять операции типа if; слыхал, что они замедляют алгоритм. Если именно в этом ошибка, то подумаю, как освобождать A и B от "балласта". (Жаль, сейчас некогда, но до ночи постараюсь успеть... а то исчезнет кнопка "Правка", и останется моя типовая ошибка тут на века... :) Если есть ещё промахи, подскажите, плз.

petr0v Тогда уж к целому положительному 64 бит 2^32 :-) да уйдёт в минус, как и положено, если к беззнаковому то переполнение разрядной сетки тут всё как всегда.

Programmist Я не советовал уйти на cuda. Просто так сказать обратил внимание. :-) Там свои колёсы, не всё паралелится, и энто изобилие процев тож ещё нужно суметь заюзать на всю катушку, эт головняк однако.

Sinus

Да именно это место. Вот сейчас и пытаюсь выпытать у Programmistа, как бесплатно по модулю 2^64 считать, без каких-либо условий.

petr0v
Где-то вы выше писали что с float медленнее получается
Медленней, конечно, если исходные данные float, но приведение их к целому, еще медленней.

что будет происходить если например к целому положительному 64 бит числу 2^64-2 прибавить 5, какой результат получится?
Число "перевернется", но за свои пределы не выйдет, такие ошибки часто бывают. Есть модули для работы с огромными числами, теоретического предела нет, только скорость сильно падает.

Sinus Да, имхо это petr0v имел ввиду, после него и я тож это просёк, а так я имел другую

S1 := SinTab(m) // Берём из таблицы два соседних значения
S2 := SinTab(m+1) // для линейной интерполяции

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

Mesh

...если к беззнаковому то переполнение разрядной сетки тут всё как всегда.

Т. е. прибавляем к целому беззнаковому 2^64-2 5 и получаем 3?

Mesh
Я не советовал уйти на cuda. Просто так сказать обратил внимание
Так поздно, теперь назад дороги нет :) В любом случае туда лезть, когда с этой прогой дело до логического конца дойдет.
Это ж как прикольно, в реалтайме на X-канальной звуковухе сразу весь эфир слушать. Жалко, что уха только два :))

petr0v Получим 2. При сложении с переполнением одна еденица из слагаемого в ноль число обратит, то есть 2^64 это 0 в беззнаковом варианте. Тут ньанс 0 занимает разряды нулями, макс число какое можно записать это 2^64-1.

Mesh

Всё же 3 получится ;)

petr0v Блин да только что проверил. Сам засумлевался. Может связано с конструкциями языка? Хотя сомнения берут, вроде всё логично. Но може где и туплю.

зы. А понял где ступил, 2^64 это уже 0, да всё верно. Я его сдуру всеми еденицами во всех битах представил. :-) Всё верно.

Mesh

Можно 4 бит сумматор взять.
2^4-1=15 максимально представимое число
2^4-2=14
14+1=15
15+1=0
0+1=1
1+1=2
2+1=3

Ну и всё вот аккумулятор готовый без лишних операций.

petr0v Ну да, тут ловушка напиши вы 0-2+5 и вопроса как бы нет. А 2^64 тут думкать надо, вот и попался. :)

Готов алгоритм без риска переполнений. Идея вот в чём. Напомню, нам нужен синус от 2пи*G*n, где n = 0,1,..., N-1, причём частота G=F/Fdiscr всегда задаётся как число, меньшее единицы. Значит, поначалу величина B=G*n заведомо меньше 1, а когда с ростом n она станет больше 1, мы из неё вычтем 1 (ибо синус от 2пи*G*n и от 2пи*(G*n-1) имеет одно и то же значение). Т.е. можем накапливать аккумулятор B, не давая ему превысить единицу. Алгоритм даже упростился:

G := F / Fdiscr // задаётся в диапазоне от 0 до 1
B := 0
Для n от 0 до N-1 с шагом 1
// Тело цикла
С := M*B
m := trunc(C) // где trunc(C) означает целую часть от С
B := B + G
Если B >= 1 то B := B - 1
S1 := SinTab^(m)
S2 := SinTab^(m+1)
OutSin⁽ⁿ⁾ := S1+(S2-S1)*(C-m)

Цикл до n=N-1 здесь только чтобы на выходе возник массив OutSin с N значениями. Если же зациклить всё безусловной петлёй, то наш алгоритм будет генерировать синус бесконечно. Как избавиться от дробных G, B, C, чтобы адрес m для элементов таблицы SinTab вычислялся только целочисленной арифметикой, пока не знаю.

Где тож косяк, ну эт уже чисто програмерский в строках с интерполяцией, так и не смог найти... У нас таблица SinTab^(m) = sin((2пи / M)*m) содержит M+1 значений и допускает m = 0,1, ... M, причём значения с m=0 и m=M одинаковы, скачка нет. А в цикле m доходит только до M-1. Имхо, всё здесь должно работать нормально. Вот, специально проверил для случая с очень грубой таблицей, М = 8:

http://s017.radikal.ru/i437/1111/c7/361252204e5a.jpg

Видим, что без интерполяции выходной массив чётко воспроизводит ступеньки из таблицы, а с интерполяцией правильно их сглаживает (здесь режим всех графиков - ступенчатый, без линейного сглаживания, но N=2048 - велико и поэтому ступенчатось собственно графиков незаметна). При G=1/4 алгоритм выдаёт OutSin = 0, 1, 0, -1, 0, ... . А при G=1/2 он выдаёт чистые нули. Если взять таблицу cos, то чётко выдаёт 1, -1, 1, -1, ... Так всё и должно быть. Если интересно, вот как наш цикл реально программируется в Маткаде - как самодельная функция, получающая извне в качестве параметров M, N, G, и массив SinTab, а возвращающая массив OutSin:

http://i073.radikal.ru/1111/65/c506d345c5e2.jpg

Всем большой спасиб. Если всё-таки чё-то не то, буду рад подсказкам. Бум учиться дальше :)

Sinus Точно, на этот ньюанс и хотел обратить внимание. Таблица должна содержать М+1 элементов. Первый и последний равны. Это нужно важный момент. В разных языках програмирования могут быть если не оговаривать это спешиал разные недоразумения. Из предыдущих описаний алгоритма это ни как не следует, и "привычное" или типовое програмирование может подложить свинью, прога будет или вылетать или генерить регулярнослучаные выборсы эт как фишка ляжет при компилировании, и как сложится что имено по этой таблице генерить будут. Ну терь когда это чётко указано, и решены проблемы с переполнениями то как и всё вроде.

А что-то в последнем релизе исчезла константа задания амплитуды? В прошлом она почему-то бралась из массива, но это несложно допереть, посему как возможную ошибку и не рассматривал.

Вобще класс, что все ошибки сами нашли.

В общем, вот такой вышел генератор:


http://www.radioscanner.ru/uploader/2011/100metrov.rar
Исходный код (в архиве) очень простой, если кто сможет оптимизировать, или добавить таблицу, буду благодарен.
Создается бинарик с разрешением 32 бита, дискрет 8000, частота 1111 Гц. Скорость работы можно измерить, но не стал добавлять ничего лишнего.
Следующим шагом, могу добавить измеритель быстродействия и ФНЧ, кстати, от ФНЧ в основной программе никакого толку нет, только ресурсы жрет, без дециматора, это пустое место.

з.ы.


А гармоники-то все равно на спектре видны, большую верность для 32-х бит, скорей всего получить не удастся :(
Вот кусочек его сигнала, конвертированного в wav: http://www.radioscanner.ru/uploader/2011/8000khzx32bit.bin.rar
С дискретом, при конвертации, ошибся в 6 раз, но это не важно.

з.ы.з.ы.
в процессорах Intel байты в многобайтных значениях переставляются так, что
младший байт идёт первым, а старший – последним
Как прочитать 4 байта Single:

function ReadSingle(b0, b1, b2, b3: Byte): Single;
var bt0, bt1, bt2, bt3 : String;
N : Longint;
S : String;
begin
bt0:= IntToHex(b3, 2); // перестановка байт и конвертация их в 2-х значный НЕХ
bt1:= IntToHex(b2, 2);
bt2:= IntToHex(b1, 2);
bt3:= IntToHex(b0, 2);
{здесь можно добавить проверку, на правильность комбинации байт}
S:= '$' + bt0 + bt1 + bt2 + bt3; // сборка числа в виде НЕХ строки ('$' в Паскале означает НЕХ)
N:=StrToInt(S); // перевод строки в целое число
Result:= Single(pointer(@N)^); // результат в виде числа с плавающей точкой, через указатель на целое число
end;

Никаких чисел с плавающей точкой нет. Все эти числа - целые числа.
При попытке конвертировать неверную комбинацию байт в число с "плавающей" точкой, возникает ошибка времени выполнения.
Подробней здесь: (http://www.delphikingdom.com/asp/viewitem.asp?catalogid=374), но это уже было.

з.ы.з.ы.з.ы.
Более "продвинутый" пример генератора с ФНЧ: http://www.radioscanner.ru/uploader/2011/100metrov2.rar
(только исходный текст).

Programmist Как вы угадали что это реал гармоники генератора а не допустим самого показометра? Просто интересно, как? Никаких чисел с плавающей точкой нет. в смысле? Вы делаете открытие очередное? По вашему ваф тег=3 ввели что б вы это опровергали? :)

Mesh
Как вы угадали что это реал гармоники генератора а не допустим самого показометра?
Это самое интересное место :) Точно смогу сказать, когда сделаю FFT, с разрешением в 128 бит.
Никаких чисел с плавающей точкой нет, в смысле?
В смысле, что вообще никаких чисел нет :) , кроме двух: "0" и "1". И всего три действия: сложение, вычитание и сдвиг, а все остальное - чудеса преобразований.
То, что мы эту точку видим на экране, как точку, вовсе не значит, что это число где-то там как дробь записано. Все хранится в целых числах. А тег=3, не знаю, зачем ввели, наверно, чтобы правильно указать читалке формат записи.

Programmist Точно смогу сказать, когда сделаю FFT, с разрешением в 128 бит. То бишь не скоро. :) Пон.

Mesh
То бишь не скоро
Кто его знает, сейчас делаю, головняк это, а иначе ответа на вопрос, откуда гармоники нет. Синус выходит 1:1, как в Куле, только быстрей немного.
Если на пальцах прикинуть, то читалка дает более 75% гармоник.

Мужики - вы уже сколько дней ударились в обсуждение цифровой обработки, воспроизведения записей...а где же обсуждение правды/лжи о влиянии "бескислородных" проводников и вообще ламповых УНЧ на качество звука? Всё ждал вердикта корифеев...а то одному тут ветерану-любителю ностальгических записей предложили кабель для динамика толщиной с палец и ценой в целую пенсию и он чуть не повёлся! Мол из особо-прокатанной меди, кристаллы чуть ли не сплошь по длине кабеля! Я как физик с трудом ему мозги подрегулировал...

думаю, уход ветки в обсуждение цифры и исчезновение в небытие адептов теплового лампового звука патефонов, разрабатывающих собственные оцифровщики пластов, подтверждают очевидное - цифра победила, и даже пласты лет 20 уже лепят с цифры, о чем все знают, и прутся от этого знания только от старых пластов

YuriVR и ничего удивительного, мне например нравиться слушать музыку с КВ с зимараниями. Вот слушать FM радио не терплю.

YuriVR и ничего удивительного, мне например нравиться слушать музыку с КВ с зимараниями. Вот слушать FM радио не терплю.

Да, читал что разработчики 174ХА42 (МС УКВ приёмника) в отличие от 174ХА34 ввели в неё генератор шума - так как при работе бесшумной настройки слушателям было "некомфортно" некоторые вообще думали что приёмник испортился и начинали его "доить" - вот и ввели ГШ в БШН для имитации "живого эфира". Психология - ничего не поделаешь!