From caf9571a694b8db80cda84f78c251c8da5ecb852 Mon Sep 17 00:00:00 2001 From: Jacob Date: Mon, 13 Feb 2023 10:35:23 +0100 Subject: [PATCH] Init --- fft_1d.py | 202 ++++++++++++++++++++++++++++++++++++++++++++++++++++ test_fft.py | 137 +++++++++++++++++++++++++++++++++++ 2 files changed, 339 insertions(+) create mode 100644 fft_1d.py create mode 100644 test_fft.py diff --git a/fft_1d.py b/fft_1d.py new file mode 100644 index 0000000..c5cf5be --- /dev/null +++ b/fft_1d.py @@ -0,0 +1,202 @@ +import numpy as np +import matplotlib.pyplot as plt +import tqdm + +import multiprocessing as mp + +RESOLUTION = 0.01 +LENGTH = 500 + + +def generate_image_from_mask(mask: np.array): + pos_mono = np.arange(0, mask.size * 2.89 * 2, 2.89) + pos_mono[::2][mask] -= 0.27 + pos_mono[1::2][mask] += 0.27 + return pos_mono + + +def pad_zero(img, length): + pad = np.zeros(length) + img = np.append(img, pad) + img = np.append(pad, img) + return img + + +def image_from_pos(pos): + length = np.max(pos) + RESOLUTION + x = np.arange(0, length, RESOLUTION) # angstrom + y = np.zeros_like(x) + ind = np.searchsorted(x, pos) + y[ind] = 1 + return y + + +def beugung(y, resolution): + fft = np.fft.fft(y) + fft_clean = np.fft.fftshift(fft) + fft_freq = np.fft.fftfreq(y.size, resolution) + fft_freq_clean = np.fft.fftshift(fft_freq) + return fft_freq_clean, np.abs(fft_clean) ** 2 + + +def gaussian_convol(img): + sigma = 100 / RESOLUTION + mu = img.size/2 + x = np.arange(0, img.size) + gauss = 1/(sigma * np.sqrt(2 * np.pi)) * \ + np.exp(- (x - mu)**2 / (2 * sigma**2)) + + return img*gauss + + +def analyisis(mask): + pos_h = generate_image_from_mask(mask) + img = image_from_pos(pos_h) + img = gaussian_convol(img) + padded = pad_zero(img, int(100 / RESOLUTION)) + freq, intens = beugung(padded, RESOLUTION) + return freq, intens + + +def get_peaks(): + orders = np.arange(1, 2, 1) + orders = orders / 5.78 + return np.array(orders) + + +def eval_peaks(freq, fft): + orders = get_peaks() + ind = np.searchsorted(freq, orders) + return fft[ind] + + +def basic_test(): + mask_h = np.zeros(LENGTH).astype(bool) + mask_l = np.ones(LENGTH).astype(bool) + + mask_mixed = np.zeros(LENGTH).astype(bool) + ind = (np.random.rand(30) * (LENGTH - 1)).astype(int) + mask_mixed[ind] = True + + mask_ner = np.zeros(LENGTH).astype(bool) + ind = (np.random.rand(1) * (LENGTH - 31)).astype(int) + ind = np.arange(ind, ind+30).astype(int) + mask_ner[ind] = True + + fig, axs = plt.subplots(4, 1) + for mask, ax in zip([mask_h, mask_l, mask_mixed, mask_ner], axs): + freq, ffty = analyisis(mask) + ax.plot(freq, ffty) + for ax in axs: + ax.plot([1.0 / 5.78, 1.0 / 2.62, 1.0 / 3.16], [0, 0, 0], "kx") + ax.plot([2.0 / 5.78, 2.0 / 2.62, 2.0 / 3.16], [0, 0, 0], "rx") + ax.plot([3.0 / 5.78, 3.0 / 2.62, 3.0 / 3.16], [0, 0, 0], "bx") + ax.set_xlim(0, 3) + plt.show() + + +def norm(arr): + return arr/np.sum(arr) + + +def norm2(arr): + return arr + # return arr/np.max(arr) + + +def next_mask(mask): + prob = np.exp((np.roll(mask, 1)*1.0 + np.roll(mask, -1)) / .1) + prob[mask] = 0.0 + prob = norm(prob) + + ind = np.random.choice(LENGTH, p=prob) + mask[ind] = True + return mask + + +def random_loop(): + mask = np.zeros(LENGTH).astype(bool) + ind = np.arange(0, LENGTH) + np.random.shuffle(ind) + percentage = [] + peaks = [] + + masks = [] + for i in ind: + mask[i] = True + freq, fft = analyisis(mask) + peak = eval_peaks(freq, fft) + percentage.append(np.mean(mask)) + peaks.append(peak) + masks.append(mask.copy()) + masks = np.array(masks) + plt.figure() + plt.imshow(masks) + plt.plot([0, 500], [406, 406]) + print() + percentage = np.array(percentage) + peaks = np.array(peaks) + return percentage, peaks + + +def nearest_loop(): + mask = np.zeros(LENGTH).astype(bool) + percentage = [] + peaks = [] + for i in range(LENGTH): + mask = next_mask(mask) + freq, fft = analyisis(mask) + peak = eval_peaks(freq, fft) + percentage.append(np.mean(mask)) + peaks.append(peak) + percentage = np.array(percentage) + peaks = np.array(peaks) + return percentage, peaks + + +def random_helper(seed): + np.random.seed(seed) + #percentage_near, peaks_near = nearest_loop() + percentage_rand, peaks_rand = random_loop() + print("done") + return percentage_rand, peaks_rand + # for i in range(peaks_near.shape[1]): + # axs[2].plot(percentage_near, norm2( + # peaks_near[:, i]), "-", label="near") + + # for i in range(peaks_rand.shape[1]): + # axs[2].plot(percentage_rand, norm2( + # peaks_rand[:, i]), ":", label="rand") + + +def random_increase(): + fig, axs = plt.subplots(3, 1) + + results = [] + for i in np.arange(10): + results.append(random_helper(i)) + + for percentage_rand, peaks_rand in results: + for i in range(peaks_rand.shape[1]): + axs[2].plot(percentage_rand, norm2( + peaks_rand[:, i]), ":", label="rand") + + for ax in [axs[0], axs[1]]: + orders = get_peaks() + ax.plot(orders, np.zeros_like(orders), "kx") + ax.set_xlim(0, 3) + + mask_l = np.ones(LENGTH).astype(bool) + mask_h = np.zeros(LENGTH).astype(bool) + freq, ffty = analyisis(mask_l) + axs[0].plot(freq, ffty) + freq, ffty = analyisis(mask_h) + axs[1].plot(freq, ffty) + plt.xlabel("percentage") + plt.ylabel("peak intensity") + plt.show() + plt.legend() + + +if __name__ == "__main__": + random_increase() diff --git a/test_fft.py b/test_fft.py new file mode 100644 index 0000000..583d6ab --- /dev/null +++ b/test_fft.py @@ -0,0 +1,137 @@ +import numpy as np +import matplotlib.pyplot as plt + + +def analysis(y, RESOLUTION): + fft = np.fft.fft(y) + fft_clean = np.fft.fftshift(fft) + fft_freq = np.fft.fftfreq(y.size, RESOLUTION) + fft_freq_clean = np.fft.fftshift(fft_freq) + + return fft_freq_clean, np.abs(fft_clean) ** 2 + + +def play_1d(): + RESOLUTION = 0.001 + LENGTH = 10000 + + x = np.arange(0, LENGTH, RESOLUTION) # angstrom + y = np.zeros_like(x) + + pos_mono = np.arange(0, x.size, 2890) + pos_mono = ( + pos_mono + np.random.normal(size=pos_mono.shape, loc=0, scale=10) + ).astype(int) + + pos_rut = np.arange(0, x.size, 5780) + + pos_rut = np.append(pos_rut, pos_rut - 3160) + + # pos_rut = (pos_rut + np.random.normal(size=pos_rut.shape, loc=0, scale=10)).astype(int) + + y[pos_rut] = 1 + y[pos_rut + 1] = 1 + y[pos_rut + 2] = 1 + + # y = np.sin(x) + + fig, axs = plt.subplots(3, 1) + ax = axs[0] + ax.plot(x, y) + + ax = axs[1] + fft_x, fft_y = analysis(y, RESOLUTION) + ax.plot(fft_x, fft_y) + ax.plot([1.0 / 5.78, 1.0 / 2.62, 1.0 / 3.16], [0, 0, 0], "kx") + ax.plot([2.0 / 5.78, 2.0 / 2.62, 2.0 / 3.16], [0, 0, 0], "rx") + ax.plot([3.0 / 5.78, 3.0 / 2.62, 3.0 / 3.16], [0, 0, 0], "bx") + ax.set_xlim(0, 3) + + +def from_mask(mask): + pos_mono = np.arange(0, mask.size * 2.89 * 2, 2.89) + + pos_mono[::2][mask] -= 0.27 + pos_mono[1::2][mask] += 0.27 + + return pos_mono + + +def image_from_pos(pos): + RESOLUTION = 0.001 + LENGTH = 1000000 + + x = np.arange(0, LENGTH, RESOLUTION) # angstrom + y = np.zeros_like(x) + ind = np.searchsorted(x, pos) + if np.any(ind > LENGTH): + print("overflow") + ind = ind[ind < LENGTH] + y[ind] = 1 + + sigma = 500 + mu = int(LENGTH / 2) + gaussian = ( + 1 + / (sigma * np.sqrt(2 * np.pi)) + * np.exp(-((x - mu) ** 2) / (2 * sigma * sigma)) + ) + # y = np.multiply(y, gaussian) + + return x, y + + +def plot_img(x, y, ax): + ax.plot(x, y) + + +if __name__ == "__main__": + RESOLUTION = 0.001 + print("Done") + + LENGTH = 1000 + mask_h = np.ones(LENGTH).astype(bool) + pos_h = from_mask(mask_h) + x, img_h = image_from_pos(pos_h) + fftx, ffty_h = analysis(img_h, RESOLUTION) + + mask_l = np.zeros(LENGTH).astype(bool) + pos_l = from_mask(mask_l) + x, img_l = image_from_pos(pos_l) + fftx, ffty_l = analysis(img_l, RESOLUTION) + print("Done") + + mask_mixed = np.zeros(LENGTH).astype(bool) + ind = (np.random.rand(400) * (LENGTH - 1)).astype(int) + mask_mixed[ind] = True + pos_mixed = from_mask(mask_mixed) + x, img_mixed = image_from_pos(pos_mixed) + fftx, ffty_mixed = analysis(img_mixed, RESOLUTION) + + print("Done") + mask_near = np.zeros(LENGTH).astype(bool) + ind = (np.random.rand(50) * (LENGTH - 1)).astype(int) + #for i in range(1, 8): + # ind = np.append(ind, ind+i) + print("Done") + + + mask_near[ind] = True + pos_near = from_mask(mask_near) + x, img_near = image_from_pos(pos_near) + fftx, ffty_near = analysis(img_near, RESOLUTION) + + fig, axs = plt.subplots(4, 1) + plot_img(fftx, ffty_h, axs[0]) + plot_img(fftx, ffty_l, axs[1]) + plot_img(fftx, ffty_mixed, axs[2]) + plot_img(fftx, ffty_near, axs[3]) + for ax in axs: + ax.plot([1.0 / 5.78, 1.0 / 2.62, 1.0 / 3.16], [0, 0, 0], "kx") + ax.plot([2.0 / 5.78, 2.0 / 2.62, 2.0 / 3.16], [0, 0, 0], "rx") + ax.plot([3.0 / 5.78, 3.0 / 2.62, 3.0 / 3.16], [0, 0, 0], "bx") + ax.set_xlim(0, 3) + # play_1d() + plt.show() + print("Done") + pass