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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.