лажаю с написанием простого DSP

[PianoIst]

Чего-то меня сегодня дернуло посмотреть ради прикола, отличается ли результат при прямой свертке в лоб и традиционный для конволюционных реверов FFT.

Открыл питон и набросал схему:

ir_l[] = массив сэмплов левого канала IR сэмпла, нормализованных 0 - 1
ir_r[] = массив сэмплов правого канала IR сэмпла, нормализованных 0 - 1

scr[] = массив моно источника

ну и дальше просто:
length = len(src) + len(ir_l)

r_l = [0] * length
r_r = [0] * length

for idx, s_frm in enumerate(src[0:100]):
    for i in range(len(ir_l)):
        v_l = ir_l[i] * s_frm
        r_l[idx + i] += v_l
        v_r = ir_r[i] * s_frm
        r_r[idx + i] += v_r
max_ = 0
for t in zip(r_l, r_r):
    for val in t:
        if max_ < abs(val):
            max_ = abs(val)

То есть, прохожусь вложенным циклом по каждому сэмплу источника, умножая его на каждый сэмпл IR и суммируя по смещению с тем, что лежит в результате.

• Получается много громче, чем оригинальный источник
• даже при нормализации (делю все это хозяйство на максимальный пик, чтобы привести в границы 0-1), получается просто шум

Правда, я не ждал, чтобы свернулся весь источник, т.к. долго. Взял только первые 100 сэмплов... Может тут причина? Но и ждать ради такого эксперимента полночи неохота)

кастую @belovw @Andrew_S.

 

[PianoIst]

свернул 10000 сэмплов) это почти весь источник

 

[Rst7]

Я вообще хз, чего там с этим питоном, как он устроен. Вот это v_l = ir_l * s_frm - это поэлементное произведение массивов? Но тогда что-то не бъется с размерностями.

В общем на псевдо-Си это выглядит так:

for(pos=ir_len; pos<(input_len-ir_len); pos++)
{
  sum=0;
  for(ipos=0; ipos<ir_len; ipos++) sum+=input[pos-ipos]*ir[ipos];
  output[pos]=sum;
}

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

 

[PianoIst]

@Rst7, спасибо! сейчас попробую.
Вот эта конструкция меня смущает

for(pos=ir_len; pos<(input_len-ir_len); pos++)
если подставляю сюда числа у меня условно получается так:
for(pos=80000; pos<(13000-80000); pos++)

и.... Или я неправильно понимаю оператор ++?

 

[Rst7]

pos++ это pos=pos+1

И да, берите размер входа больше длины импульса, а то фигня получится. Ну либо вход нулями добивайте сначала.
[DOUBLEPOST=1554244091][/DOUBLEPOST]

for(pos=80000; pos<(13000-80000); pos++)

Такой цикл ни разу не выполнится. Потому что 80000 (длина импульса) превышает размер входа (13000). Так что надо на входе сначала положить 80000 нулей, а потом уже интересующий сигнал, тогда правильно будет сворачивать.
[DOUBLEPOST=1554244148][/DOUBLEPOST]Повторюсь, свертка - она в обратную сторону, потому чтобы самый первый семпл на входе обработать - надо знать все предыдущие. Ну обычно считают, что они все нули.
[DOUBLEPOST=1554244413][/DOUBLEPOST]Вариант без добивания нулями

for(pos=0; pos<input_len; pos++)
{
  sum=0;
  for(ipos=0; ipos<ir_len; ipos++) if (ipos<=pos) sum+=input[pos-ipos]*ir[ipos];
  output[pos]=sum;
}

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

 

[PianoIst]

@Rst7, простите, чет я не догоняю...
Прямой перенос дает следующий результат

Я даже себе картинку нарисовал как понимаю должно происходить, вроде все сходится, кроме длины первого цикла:

Но сумма достаточно резко улетает в облака, и почти не прогибается вниз по дороге. Вот сейчас у меня уже сворачивается 54900 сэмпл, где инпута в помине нет (нулями в конце добил для этой итерации), а сумма все растет
может надо среднее арифметическое брать?)

 

[Rst7]

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

А посмотрите, у Вас там случайно при чтении файлов не получилось unsigned int вместо signed int? Как файлы читаете и преобразуете во float?

 

[PianoIst]

не получилось unsigned int вместо signed int?

именно.
Почему-то когда пытаюсь приводить из unsigned int через (sample - 32768) / 32768 и соответственно обратно через sample * 32768 + 32768 получается повышение уровня в 2 раза...

то есть сейчас я wav распаковываю как unsigned и упаковываю как unsigned, и работаю тоже в диапазоне 0 - 1.
Самое дурацкое, что когда я просто пытаюсь диапазон сместить на 0.5 и умножить все на 2 тоже лажа выходит какая-то...
Есть подозрение, что у меня где-то сэмплы идут как ссылки, а не как данные, но... хз где)

 

[Rst7]

Ну это же неправильно. wav - это signed, и преобразовывать надо в диапазон -1...1. Не забывая перед преобразованием при записи проверить на диапазон -1<=val<=1 и ограничить до умножения на 32767.

 

[H-ron]

@PianoIst, точный листинг скиньте, плз.

 

[Rst7]

float val;

int inputsample;

int outputsample;

туда: 

val=inputsample;
val/=32767.0;

оттуда:

if (val>1.0) val=1.0;
if (val<-1.0) val=-1.0;
val*=32767.0;
outsample=val;

 

[PianoIst]

@Rst7, я балда... Скопировал пример с библиотекой распаковки байтов, и не посмотрел, что там аргумент стоит uint16...
Вечером перепишу, отпишусь, должно заработать)
@Hron, наверное, сейчас пока не актуально (у меня скоро урок, потом поправлю, проверю)
но на всякий:

Спойлер: код

import wave
import numpy
import struct
range_ = 32768 * 2

with wave.open('IR.wav', 'r') as ir:
    print(ir.getparams())
    source = ir.readframes(ir.getnframes())
    ir_frames = [
        (frm / range_)
        for frm in numpy.frombuffer(source,
                                    dtype='uint16')
    ]
    print(len(ir_frames), len(ir_frames) / 2)
    ir_l = list()
    ir_r = list()
    for i in range(int(len(ir_frames) / 2)):
        ir_l.append(ir_frames[i * 2])
        ir_r.append(ir_frames[i*2 + 1])

with wave.open('source16.wav', 'r') as src:
    print(src.getparams())
    source = src.readframes(src.getnframes())
    src_frames = [
        (frm / range_)
        for frm in numpy.frombuffer(source,
                                    dtype='uint16')
    ]
# print([frm for frm in numpy.frombuffer(source, dtype='uint16')])
# print(src_frames)

# src_frames = [0] * len(ir_l) + src_frames
src_frames += [0] * len(ir_l)
length = len(src_frames) + len(ir_l)
print(len(src_frames), len(ir_l), len(ir_r), length)

r_l = [0] * length
r_r = [0] * length

print(
    'first: ',
    len(ir_l),
    len(src_frames) - len(ir_l),
    len(src_frames[:100])
)
for pos in range(len(src_frames)):
    # pos = c_pos + len(ir_l)
    sum_ = 0
    if pos % 100 == 0:
        print(pos, sum_)
    for ipos in range(len(ir_l)):
        if ipos > pos:
            continue
        sum_ += src_frames[pos - ipos] * ir_l[ipos]
        # print('   ', src_frames[pos - ipos], ir_l[ipos], sum_)
    r_l[pos] = sum_
    if pos % 100 == 0:
        print(pos, sum_)

# breakpoint()

max_ = 1
for t in zip(r_l, r_r):
    for val in t:
        if max_ < abs(val):
            max_ = abs(val)
# print(len(r_l), r_l[88248])
max_ += 0.1
print(max_)

with wave.open('out.wav', 'w') as out:
    out.setnchannels(1)
    out.setnframes(length + 1)
    out.setsampwidth(2)
    out.setframerate(44100)
    # print((src_frames[1] + 1) * (32767/2))
    for frm_l, frm_r in zip(r_l, r_r):
        val_l = int(frm_l * range_ / max_)
        val_r = int(frm_r * range_ / max_)
        # if val_l > 0xffff:
        #     val_l = 0xffff
        # if val_r > 0xffff:
        #     val_r = 0xffff
        print(out.tell(), val_l, val_r)
        out.writeframesraw(struct.pack('H', val_l))
        # out.writeframesraw(struct.pack('H', val_r))