Контакты

Вы можете связаться со мной, если это необходимо


Семестр VII. Рентгеноструктурный анализ (часть 3)

Восстановление функции по коэффициентам ее ряда Фурье

Задача

Исследовать качество восстановления функции одной переменной в зависимости от того, какие и сколько гармоник ряда Фурье используются для ее восстановления.
Для этого необходимо симулировать эксперимент по получению данных об амплитудах и фазах гармоник Фурье и восстановить по ним функцию ЭП.

Модель для компьютерного эксперимента

  • Для симуляции использовалась модель одномерной электронной плотности.
  • На отрезке от 0 до 30 (ангстрем) расположены две молекулы. Атомы в молекуле связаны ковалентно и находятся на расстоянии 1-1.5 анстрем друг от друга. Молекулы расположены на расстоянии 3-5 ангстрем (водородная связь или гидрофобное взаимодействие между ними).
  • Функция имеет вид нескольких гауссовых кривых с центром в разных точках. Электронные плотности (ЭП) атомов описываются этими гауссовыми кривыми, которые представляют собой пики на графике.
  • Гауссова функция определяется числами lambda,beta,gamma (параметрами) по формуле: gauss = lambda*exp(-(beta^2)*(X-gamma)^2)
  • Параметры гауссовских кривых соответствуют различным атомам. Различная высота пиков соответствует различному количеству электронов в атоме, ширина пика примерно соответствует диаметру атома, то есть порядка 1 ангстрема(электронная плотность резко спадает на таком расстоянии от центра).

Симуляция

Функция электронной плотности атомов на этом отрезке была задана с помощью скрипта compile-func.py (из пяти пиков - 3+2, всего две "молекулы" - модель водородной связи (3.5) ):
python compile-func.py -g 16,3.5,4.72+12,3.5,6.54+14,3.5,7.82+1,3,11.32+14,3.5,12.33
Каждое "слагаемое" в формуле для параметра -g соответствует одному атому (пику) и содержит в себе параметры, разделенные запятыми:
  • высота пика (число электронов)
  • ширинa пика (размер атома)
  • положение пика на прямой (так как симуляция одномерная)
Выходной файл скрипта compile-func.py - текстовый файл func.txt, в котором описана заданная функция (в виде пар X – Y).
Рис.1. Функция электронной плотности атомов двух воображаемых молекул на отрезке [0;30]
Экспериментальной информацией, на самом деле, являются данные об амплитудах и фазах гармоник Фурье. Коэффициенты Фурье можно рассчитать при помощи скрипта func2fourier.py. В результате получается файл вида garm.txt:
#no.	F	phi
0	4.72654359376	0
1	1.39720718995	-2.12035906411
2	1.98245142474	-1.43298013822
3	0.433019796648	1.19952810495
4	1.16454822148	-2.28301295427
5	0.827520704547	-2.09439511117

python func2fourier.py -i func.txt -o garm.txt
Выходной файл имеет формат:
<номер гармоники>  <амплитуда>  <фаза>
Добавление гауссового шума к амплитудам (параметр -F <число>) и фазам (-P <число>) искажает все вычисленные амплитуды и фазы. Пример: -F 20 (шум 20%) приводит к тому, что к каждой амплитуде прибавляется случайное число, распределенное нормально с параметрами: среднее = 0, среднее квадратичное отклонение (сигма)=0.2*F. Аналогично действует параметр -P <число>. Это позволяет еще более точно приблизится к эксперименту.

Восстановление функции электронной плотности. Полные наборы гармоник.

Для сравнения восстановленной функции с исходной использовались следующие категории:
  • Отличное восстановление – по графику восстановленной функции можно определить положение максимума всех гауссовых слагаемых функции ("атомов").
  • Хорошее восстановление – можно угадать положение всех максимумов, зная число слагаемых ("атомов"), хотя на восстановленной функции максимумы от атомов не отличимы от шума.
  • Среднее восстановление – положение каких-то атомов определить по восстановленной функции нельзя, других - можно.
  • Плохое восстановление – положение атомов определить не представляется возможным; можно только предсказать примерный размер "молекулы".
Коэффициенты ряда Фурье рассчитаны при помощи скрипта func2fourier.py. Результат - файл garm_full.txt:
python func2fourier.py -i func.txt -o garm_full.txt -F 0 -P 0
Из 499 гармоник отбиралось n0 гармоник (полный набор - 0, 1, ..., n0), для которых по графику восстановленной по ним функции можно было определить положение максимумов всех гауссовых слагаемых («отличное восстановление»):
python fourier-filter.py -r 0-15 -i garm_full.txt -o garm_full_0-15.txt
python fourier2func.py -f func.txt -i garm_full_0-15.txt -o recr_func_0-15.txt -s
Исходная функция n0 = 5 (плохое) n0 = 10 (плохое)
Исходная функция n0 = 15 (плохое) n0 = 20 (среднее)
Исходная функция n0 = 25 (среднее) n0 = 30 (хорошее)
Исходная функция n0 = 35 (хорошее) n0 = 40 (хорошее)
Исходная функция n0 = 45 (отличное) n0 = 70 (отличное)


Таким образом, четыре пика для тяжелых атомов хорошо различимы уже на n = 20, но пик легкого атома выделяется только при n = 45, поэтому за минимальный полный набор гармоник примем n0=45. Таким образом, основная трудность - разрешение соседних пиков, а также отличие пиков малых атомом (водород) от шума.

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

Добавление шума

Теперь добавим шум к амплитудам и фазам. По полному набору гармоник (0-45) были построены графики восстановленной функции при ненулевом значении шума.
# Добавление шума к амплитуде
python func2fourier.py -F 20 -i func.txt -o ampl_noise20.txt
python fourier-filter.py -r 0-45 -i ampl_noise20.txt -o ampl_noise20_45.txt
python fourier2func.py -i ampl_noise20_45.txt -o recr_func_ampl_noise20_45.txt -s

# Добавление шума только к фазе
python func2fourier.py -P 20 -i func.txt -o phase_noise20.txt
python fourier-filter.py -r 0-45 -i phase_noise20.txt -o phase_noise20_45.txt
python fourier2func.py -i phase_noise20_45.txt -o recr_func_phase_noise20_45.txt -s

# Добавление шума к амплитуде и к фазе
python func2fourier.py -F 20 -P 20 -i func.txt -o both_noise20.txt
python fourier-filter.py -r 0-45 -i both_noise20.txt -o both_noise20_45.txt
python fourier2func.py -i both_noise20_45.txt -o recr_func_both_noise20_45.txt -s
Исходная функция Нет шума Шум амплитуды 20 %
Исходная функция Шум фазы 20 % Шум амплитуды и фазы 20 %


Из рисунков видно, что ошибка в фазе больше "мешает" верно определить пики, чем такой же шум в амплитуде. Поэтому решение проблемы фаз в настоящем эксперименте - важная задача, так как неверное определение фазы сильно влияет на результат.

Неполные наборы гармоник. Удаление первых гармоник

Удаление первых (0, 0-1, 0-3) гармоник поможет понять, как первые гармоники влияют на качество результата.
# Удаление первой гармоники
python fourier-filter.py -r 1-45 -i garm_full.txt -o garm_1-45.txt
python fourier2func.py -i garm_1-45.txt -o recr_func_1-45.txt -s

# Удаление первых двух гармоник
python fourier-filter.py -r 2-45 -i garm_full.txt -o garm_2-45.txt
python fourier2func.py -i garm_2-45.txt -o recr_func_2-45.txt -s

# Удаление первых трех гармоник
python fourier-filter.py -r 3-45 -i garm_full.txt -o garm_3-45.txt
python fourier2func.py -i garm_3-45.txt -o recr_func_3-45.txt -s

Исходная функция Гармоники 0-45 Гармоники 1-45
Исходная функция Гармоники 2-45 Гармоники 3-45


При удалении первой гармоники функция смещается по оси ординат; при этом результат восстановления можно считать «отличным». При удалении первых двух и трех гармоник восстановленная функция сильно искажается, но значения максимумов для больших атомов можно однозначно определить (водород, как и ожидается, не отличим от шума). Такомобразом, удаление первых гармоник значительно искажает фоновое значение (фон теперь не образует "прямую" вдоль оси абсцисс) и увеличивается шум.

Неполные наборы гармоник. Удаление гармоник из середины набора

Было решено удалить 6,7% и 11,1% гармоник из середины набора (гармоники 16, 24 и 32 (6%), 8, 40 (12%)). Такое удаление гармоник не дает никакого радикального эффекта, хоть разрешение, а значит и качество восстановления, снижается вследствие неполноты данных (можно охарактеризовать как «среднее»). Таким образом, удаление гармоник в середине набора приводит к тому,что максимумы от «атомов» становятся менее отличимы от фонового шума, вследствие его увеличения.
# Удаление 6,7% гармоник
python fourier-filter.py -r 0-15,17-23,25-31,33-45 -i garm_full.txt -o garm_6.txt
python fourier2func.py -i garm_6.txt -o recr_func_6.txt -s

# Удаление 11,1% гармоник
python fourier-filter.py -r  0-7,9-15,17-23,25-31,33-39,41-45 -i garm_full.txt -o garm_11.txt
python fourier2func.py -i garm_11.txt -o recr_func_11.txt -s
Исходная функция Полный набор гармоник (0...45)
6,7% гармоник удалено 11,1% гармоник удалено

Добавление гармоники n0+10

К набору гармоник 0-45 добавлена гармоника 55, что ни в коей мере не изменило качество восстановления.
# Добавление 55 гармоники
python fourier-filter.py -r 0-45,55 -i garm_full.txt -o garm_55.txt
python fourier2func.py -i garm_55.txt -o recr_func_55.txt -s

Исходная функция Полный набор гармоник (0...45) Гармоники (0...45, 55)



Таблица 1. Восстановление функции по коэффициентам ряда Фурье

Набор гармоник Разрешение, A Полнота данных,% Шум амплитуды (% от F) Шум фазы (% от phi) Качество восстановления
Полный набор гармоник
0-5 6 100 0 0 Плохое
0-10 3 100 0 0 Плохое
0-15 2 100 0 0 Плохое
0-20 1.5 100 0 0 Среднее
0-25 1.2 100 0 0 Среднее
0-30 1 100 0 0 Хорошее
0-35 0.86 100 0 0 Хорошее
0-40 0.75 100 0 0 Хорошее
0-45 0.67 100 0 0 Отличное
0-45 0.67 100 20 0 Хорошее
0-45 0.67 100 0 20 Среднее
0-45 0.67 100 20 20 Среднее
Неполный набор гармоник
1-45 0.67 97.82 0 0 Отличное
2-45 0.67 95.65 0 0 Хорошее
3-45 0.67 93.48 0 0 Хорошее
0-15,17-23,25-31,33-45 0.67 93.48 0 0 Среднее
0-7,9-15,17-23,25-31,33-39,41-45 0.67 89.13 0 0 Среднее
0-45, 55 0.67 100 0 0 Отличное