Что я делаю неправильно при преобразовании этого кода MATLAB в Python?

Question

Я работаю над преобразованием этого кода MATLAB, который генерирует форму волны, в Python. Для контекста это моделирование реакции возбуждения полосы с помощью силового микроскопа atomi c (не имеет отношения к ошибке кода). График, который создается из r_ve c в MATLAB, отличается от того, когда я генерирую его в Python, правильно ли я конвертирую код MATLAB в Python?

Я прикрепил график r_vec это построено в MATLAB, и именно так оно должно выглядеть по сравнению с тем, что генерирует Python. Вот код MATLAB ниже. После этого мой преобразованный код Python. Я подозреваю, что проблема связана с для l oop, поскольку все переменные, которые я подтвердил, равны между кодом MATLAB и Python, хотя я не могу быть уверен, и не я был тем, кто создал Код MATLAB. Я потратил часы, пытаясь выяснить, в чем может быть проблема, но я, кажется, нашел ее.

График MATLAB

Python Plot

Python участок

Код MATLAB

%simulate BE response over a line scan

% define experimental parameters
IO_rate = 4E6; %[samples/sec]
N_pixels = 128; % number of pixels along a line scan
N_points_per_pixel = 2^13; % number of data points per pixel
w1 = 200E3; % lower edge of band
w2 = 400E3; % upper edge of band
noise_level = .1; %add noise to the signal

w_vec = -IO_rate/2: IO_rate/N_points_per_pixel : IO_rate/2-IO_rate/N_points_per_pixel;

p_vec = (0:N_pixels-1)/N_pixels;
A_vec = sin(2*pi*3*p_vec)+2; %amplitude
wo_vec = 250E3 + 100E3*p_vec; %resonance
Q_vec = 100 - 50*p_vec; % Q-factor
phi_vec = sign(p_vec-.5); % phase
D_vec = zeros(size(w_vec));

D_vec( ((abs(w_vec)<w2) + (abs(w_vec)>w1)) == 2 ) = 1; % drive bins located within upper and lower band edges
band_ind = find( (((w_vec)<w2) + ((w_vec)>w1)) == 2 );

d_vec = fftshift(ifft(ifftshift(D_vec))); % find drive signal in the time domain

R_mat = zeros(N_pixels,N_points_per_pixel);
r_mat = zeros(N_pixels,N_points_per_pixel);

for k1 = 1 : N_pixels
    H_vec = (A_vec(k1).*wo_vec(k1).^2).*exp(1i*phi_vec(k1))./(w_vec.^2 + 1i*wo_vec(k1)*w_vec/Q_vec(k1) - wo_vec(k1).^2); %cantilever transfer function
    R_mat(k1,:) = (H_vec.*D_vec); %response of the cantilever in the Fourier domain
    r_mat(k1,:) = fliplr((real((ifft(fftshift(R_mat(k1,:)))))));
end

r_vec = reshape(r_mat.',[ 1 N_pixels*N_points_per_pixel]);

Python Код

IO_rate = 4000000 #[samples/sec]
N_pixels = 128 # number of pixels along a line scan
N_points_per_pixel = 8192 # number of data points per pixel
w1 = 200000 # lower edge of band
w2 = 400000 # upper edge of band
noise_level = .1 #add noise to the signal

w_vec = np.arange(-IO_rate/2, IO_rate/2-IO_rate/N_points_per_pixel + 1, IO_rate/N_points_per_pixel) #frequency vector over a pixel

p_vec = np.arange(0, N_pixels+1)/N_pixels
p_vec = np.delete(p_vec, -1) #delete last element
A_vec = np.sin(2*np.pi*3*p_vec)+2 #amplitude
wo_vec = 250000 + 100000*p_vec #resonance
Q_vec = 100 - 50*p_vec # Q-factor
phi_vec = np.sign(p_vec-.5) # phase

D_vec = np.zeros((np.size(w_vec)), dtype=complex)
D_vec = np.real(D_vec)

ind = (abs(w_vec)<w2) & (abs(w_vec)>w1)
D_vec[ind] = 1 #assign those indices to 1.
band_ind = np.argwhere(((w_vec)<w2) & ((w_vec)>w1)) #check again
band_ind = np.delete(band_ind, [0]) #delete first element to correspond exactly w/ MATLAB
band_ind = np.append(band_ind, 4916) #add 4916 as last element to correspond w/ MATLAB

d_vec = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(D_vec)))
d_vec = np.real(d_vec)

R_mat = np.zeros((N_pixels,N_points_per_pixel))
r_mat = np.zeros((N_pixels,N_points_per_pixel))

for k1 in range(N_pixels):
    H_vec = ((A_vec[k1]*wo_vec[k1]**2)*np.exp(1j*phi_vec[k1])/(w_vec**2 + 1j*wo_vec[k1]*w_vec/Q_vec[k1] - wo_vec[k1]**2)) #cantilever transfer function
    R_mat[k1,:] = (H_vec*D_vec) #re  sponse of the cantilever in the Fourier domain
    r_mat[k1,:] = np.flip((np.real((np.fft.ifft(np.fft.fftshift(R_mat[k1,:]))))))

r_vec = np.reshape(r_mat.transpose(), (N_pixels*N_points_per_pixel))

Answer 1

Вы можете переписать свой код Python, используя трансляцию

import matplotlib.pyplot as plt
from numpy import (abs, arange, exp, linspace, logical_and, pi,
                   sign, sin, where)
from numpy.fft import fft, fftshift, ifft, ifftshift

rate = 2*2E6
npx = 128
nptpx = 8192
ω_low, ω_hi = 200000, 2*200000

ω, step = linspace(-rate/2, rate/2, nptpx,
                       endpoint=0, retstep=1)
p = arange(npx)/npx
A = sin(6*pi*p)+2
ω0 = 250000+100000*p
Q = 100-50*p
φ = sign(p-1/2)
D = where(logical_and(abs(ω)>ω_low, abs(ω)<ω_hi), 1, 0)
R = (D * (A * ω0**2 * exp(1j*φ))[:,None] /
     (ω**2 + 1j*ω*(ω0/Q)[:,None] - (ω0**2)[:,None]))

r = ifft(fftshift(R, axes=1))[:,::-1].real.reshape(npx*nptpx)

plt.figure(figsize=(12,4), constrained_layout=1)
plt.plot(r, lw=.4)
plt.xlim(0,len(r)-1)
plt.show()

Compare with your plot

График MATLAB