Детерминированная деконволюция Фурье
Я не думаю, что в Scala есть что-то особенное или сложное в фигурных скобках. Чтобы справиться с кажущимся сложным использованием их в Scala, просто сохраните пару простых вещей:
- фигурные скобки образуют блок кода, который оценивает последнюю строку кода (почти все языки делают это)
- функция, если требуется, может быть сгенерирована блоком кода (следует правилу 1)
- фигурные скобки могут быть опущены для однострочного кода, за исключением случая (выбор Scala)
- круглые скобки могут быть опущены при вызове функции с кодовым блоком в качестве параметра (выбор Scala)
Давайте объясним несколько примеров в трех предыдущих правила:
val tupleList = List[(String, String)]()
// doesn't compile, violates case clause requirement
val filtered = tupleList.takeWhile( case (s1, s2) => s1 == s2 )
// block of code as a partial function and parentheses omission,
// i.e. tupleList.takeWhile({ case (s1, s2) => s1 == s2 })
val filtered = tupleList.takeWhile{ case (s1, s2) => s1 == s2 }
// curly braces omission, i.e. List(1, 2, 3).reduceLeft({_+_})
List(1, 2, 3).reduceLeft(_+_)
// parentheses omission, i.e. List(1, 2, 3).reduceLeft({_+_})
List(1, 2, 3).reduceLeft{_+_}
// not both though it compiles, because meaning totally changes due to precedence
List(1, 2, 3).reduceLeft _+_ // res1: String => String = <function1>
// curly braces omission, i.e. List(1, 2, 3).foldLeft(0)({_ + _})
List(1, 2, 3).foldLeft(0)(_ + _)
// parentheses omission, i.e. List(1, 2, 3).foldLeft(0)({_ + _})
List(1, 2, 3).foldLeft(0){_ + _}
// block of code and parentheses omission
List(1, 2, 3).foldLeft {0} {_ + _}
// not both though it compiles, because meaning totally changes due to precedence
List(1, 2, 3).foldLeft(0) _ + _
// error: ';' expected but integer literal found.
List(1, 2, 3).foldLeft 0 (_ + _)
def foo(f: Int => Unit) = { println("Entering foo"); f(4) }
// block of code that just evaluates to a value of a function, and parentheses omission
// i.e. foo({ println("Hey"); x => println(x) })
foo { println("Hey"); x => println(x) }
// parentheses omission, i.e. f({x})
def f(x: Int): Int = f {x}
// error: missing arguments for method f
def f(x: Int): Int = f x
1 ответ
1) Как вы правильно поняли, требование к масштабированию связано с дискретным преобразованием Фурье. Лучший способ получить это - вычислить деконволюцию двух унифицированных единичных сигналов. Их DFT равен n 0 0 0 ...., где n - количество точек DFT. Следовательно, отношение r_f равно 1 0 0 0 0, и его обратное БПФ, вычисленное с помощью np.fft.ifft(), равно 1 / n 1 / n 1 / n ... Правильный сигнал, полученный в результате деконволюции, должен был быть 1 / T 1 / T 1 / Т ..., где Т = 10. это длина кадра.
Как следствие, правильное масштабирование для выполнения деконволюции составляет n/T= len(r_f)/10.
r_if=r_if*len(r_if)/10.
2) Деконволюционный сигнал транслируется на половину периода. Это связано с тем, что ядро Гаусса центрируется по центру кадра. Просто сдвиньте ядро на половину периода, и проблема решена. Для этого можно применить функцию np.fft.fftshift():
f_yb = np.fft.fft(np.fft.fftshift(y_b)) #fft the filter
РЕДАКТИРОВАТЬ: Чтобы исследовать причину перевода, давайте сосредоточимся на случае, когда ядро деконволюции является очень узким гауссовым распределением, почти соответствующим Распределение Дирака. Ваш входной сигнал представляет собой гауссову кривую с центром в нуле, а кадр дискретизируется между -5 и 5. Аналогично, ядро деконволюции представляет собой дираковский центр с центром в нуле. Как следствие, деконволюция сигнала должна быть идентична входному сигналу: гауссова кривая с центром в нуле. Тем не менее, DFT, реализованный в FFTW и, следовательно, np.fft.fft() , вычисляется как кадр, начинающийся с 0 и заканчивающийся на 10 , который выбирается в точках 10j / n, где j находится в [0..n- 1], частоты в пространстве Фурье равны k / 10, где k в [0..n / 2, -n / 2 + 1 ..- 1]. Как следствие, этот ДПФ воспринимает ваш сигнал как гауссиан с центром в 5, а ядро деконволюции - как дираковский с центром в 5. Свертка функции f (t) с дельтой Дирака (t-t0) с центром в точке t0 является просто переведенная функция f (t-t0). Следовательно, результатом деконволюции, вычисленной в np.fft.fft(), является входной сигнал, переведенный на половину периода. Поскольку входной сигнал центрируется на 0 в кадре [-5,5], выходной сигнал, вычисленный с помощью np.fft.fft(), центрируется на -5 (или эквивалентно 5 из-за периодичности). Сдвиг ядра устраняет несоответствие между тем, что мы думаем о кадре как [-5 5], и np.fft.ifft() обрабатываем его так, как если бы это было [0 10].
Фильтры часто предназначены для уменьшения влияния высокочастотных шумов. Таким образом, деконволюция вызывает потенциальное увеличение высокочастотного шума. Скрининг частот, как вы сделали, является потенциальным решением. Обратите внимание, что это точно эквивалентно свертке сигнала с определенным фильтром!
В диапазоне томографической реконструкции алгоритм отфильтрованного обратного проецирования требует применения линейного фильтра, который значительно раздувает высокочастотный шум. Здесь предложен фильтр Винера : этот тип фильтра может быть спроектирован так, чтобы минимизировать среднеквадратичную ошибку на деконволюционном сигнале, учитывая SNR извилистого сигнала. Тем не менее, требуется некоторое предположение относительно спектральной плотности мощности сигнала и шума.