1

I am working to convert this MATLAB code that generates a waveform to Python. For context this is a simulation of band excitation response from an atomic force microscope (not relevant to code error). Here is the MATLAB code

%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; %frequency vector over a pixel

% vary A, wo, Q, and phase over pixels
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

% build drive signal, define in the Fourier domain
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

% build response at each pixel
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
    
    %determine response in the time domain (this is a little hokey, but it should work for simulation)    
    r_mat(k1,:) = fliplr((real((ifft(fftshift(R_mat(k1,:)))))));    
end

% build full response in the time domain;
r_vec = reshape(r_mat.',[ 1 N_pixels*N_points_per_pixel]);

% add noise
r_vec = r_vec + noise_level*2*(rand(size(r_vec))-.5);

%save response as a .mat (which can be read into python if needed)

Here is what I have so far for converting this to python code

#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 = 8192; # 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 = np.arange(-IO_rate/2, IO_rate/2-IO_rate/N_points_per_pixel + 1, IO_rate/N_points_per_pixel)
# vary A, wo, Q, and phase over pixels
p_vec = np.arange(0, N_pixels-1)/N_pixels
A_vec = np.sin(2*np.pi*3*p_vec)+2 #amplitude
wo_vec = 250E3 + 100E3*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))
ind = (abs(w_vec)<w2) & (abs(w_vec)>w1);
D_vec[ind] = 1; #assign those indices to 1.
band_ind = np.nonzero(((w_vec)<w2) & ((w_vec)>w1));

d_vec = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(D_vec))) #find drive signal in the time domain
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

After executing what I have so far in the for loop I get TypeError: 'numpy.ndarray' object is not callable so I'm not sure what I'm doing wrong?

1
  • 2
    Do you understand the difference between MATLAB and python when it comes to indexing? The use of () versus []? Do you know what a traceback is? And how to identify which line and variable is the problem? Commented Jul 23, 2020 at 5:20

2 Answers 2

1

The issue exists in the indexing of the vector in the loop.
The code should be :

    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

Also there seems to be an issue in the loop. Did you mean to write:

for k1 in range(N_pixels-1):
      
Sign up to request clarification or add additional context in comments.

3 Comments

Why would you suggest OP loop over all elements except the first and last? The MATLAB code loops over the full array.
@CrisLuengo How could I write it so that it loops the full array instead?
@RyanF: your for k1 in range(N_pixels): is correct.
0
  1. The reason is that you use the () operator instead of [] to access the items in an array (i.e. you use the MatLab style of element referencing instead of a Phytonian style).
  2. Also, in the np.arrange() you should pass N_pixels, otherwise you'll get an IndexError: index 127 is out of bounds for axis 0 with size 127 error.
  3. As Python uses indentation as separation between lines, there's no need to add ; at the end of each line.

Here is the corrected version:

import numpy as np
#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 = 8192 # 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 = np.arange(-IO_rate/2, IO_rate/2-IO_rate/N_points_per_pixel + 1, 
IO_rate/N_points_per_pixel)
# vary A, wo, Q, and phase over pixels
p_vec = np.arange(N_pixels)/N_pixels
A_vec = np.sin(2*np.pi*3*p_vec)+2 #amplitude
wo_vec = 250E3 + 100E3*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))
ind = (abs(w_vec)<w2) & (abs(w_vec)>w1)
D_vec[ind] = 1; #assign those indices to 1.
band_ind = np.nonzero(((w_vec)<w2) & ((w_vec)>w1))

d_vec = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(D_vec))) #find drive signal in the time domain
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

Cheers.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.