0

I have successfully converted a MATLAB source code to Python - but the plot output are just not matching. I have doubly verified the values of each variables bot in Python and Octave - both of them are same also.

Octave Plot Output:

enter image description here

Python Matplotlib Output:

enter image description here

Octave Code:

clear
N  = 10^3; % number of symbols
am = 2*(rand(1,N)>0.5)-1 + j*(2*(rand(1,N)>0.5)-1); % generating random binary sequence
fs = 10; % sampling frequency in Hz

% defining the sinc filter
sincNum = sin(pi*[-fs:1/fs:fs]); % numerator of the sinc function
sincDen = (pi*[-fs:1/fs:fs]); % denominator of the sinc function
sincDenZero = find(abs(sincDen) < 10^-10);
sincOp = sincNum./sincDen;
sincOp(sincDenZero) = 1; % sin(pix/(pix) =1 for x =0

% raised cosine filter
alpha = 0.5;
cosNum = cos(alpha*pi*[-fs:1/fs:fs]);
cosDen = (1-(2*alpha*[-fs:1/fs:fs]).^2);
cosDenZero = find(abs(cosDen)<10^-10);
cosOp = cosNum./cosDen;
cosOp(cosDenZero) = pi/4;

gt_alpha5 = sincOp.*cosOp;

alpha = 1;
cosNum = cos(alpha*pi*[-fs:1/fs:fs]);
cosDen = (1-(2*alpha*[-fs:1/fs:fs]).^2);
cosDenZero = find(abs(cosDen)<10^-10);
cosOp = cosNum./cosDen;
cosOp(cosDenZero) = pi/4;
gt_alpha1 = sincOp.*cosOp;

% upsampling the transmit sequence
amUpSampled = [am;zeros(fs-1,length(am))];
amU = amUpSampled(:).';

% filtered sequence
st_alpha5 = conv(amU,gt_alpha5);
st_alpha1 = conv(amU,gt_alpha1);

% taking only the first 10000 samples
st_alpha5 = st_alpha5([1:10000]);
st_alpha1 = st_alpha1([1:10000]);

st_alpha5_reshape = reshape(st_alpha5,fs*2,N*fs/20).';
st_alpha1_reshape = reshape(st_alpha1,fs*2,N*fs/20).';

close all
figure;
st_alpha5_reshape
plot([0:1/fs:1.99],real(st_alpha5_reshape).','b');  
title('eye diagram with alpha=0.5');
xlabel('time')
ylabel('amplitude')
axis([0 2 -1.5 1.5])
grid on

figure;
plot([0:1/fs:1.99],real(st_alpha1_reshape).','b'); 
title('eye diagram with alpha=1')
xlabel('time')
ylabel('amplitude')
axis([0 2 -1.5 1.5 ])
grid on

Python code:

j = np.complex(0,1)
N = 10**3
#% number of symbols
am = 2.*(np.random.rand(1., N) > 0.5)-1.+np.dot(j, 2.*(np.random.rand(1., N) > 0.5)-1.)
#% generating random binary sequence
fs = 10.
#% sampling frequency in Hz
#% defining the sinc filter
sincNum = np.sin(np.dot(np.pi, np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs))))))
#% numerator of the sinc function
sincDen = np.dot(np.pi, np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs)))))
#% denominator of the sinc function
sincDenZero = np.where(abs(sincDen) < 10**-10);
sincOp=sincNum/sincDen
sincOp[int(sincDenZero[0])-1] = 1.
#% raised cosine filter
alpha = 0.5
cosNum = np.cos(np.dot(np.dot(alpha, np.pi), np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs))))))
cosDen = 1.-np.dot(2.*alpha, np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs)))))**2.
cosDenZero = np.where(abs(cosDen)<10**-10);
cosOp = cosNum/cosDen
cosOp[int(cosDenZero[0][0])-1] = np.pi/4.
cosOp[int(cosDenZero[0][1])-1] = np.pi/4.
gt_alpha5 = sincOp*cosOp
alpha = 1.
cosNum = np.cos(np.dot(np.dot(alpha, np.pi), np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs))))))
cosDen = 1.-np.dot(2.*alpha, np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs)))))**2.
cosDenZero = np.where(abs(cosDen)<10**-10);
cosOp = cosNum/cosDen
cosOp[int(cosDenZero[0][0])-1] = np.pi/4.
cosOp[int(cosDenZero[0][1])-1] = np.pi/4.
gt_alpha1 = sincOp*cosOp
#% upsampling the transmit sequence
#amUpSampled = np.array(np.vstack((np.hstack((am)), np.hstack((np.zeros((fs-1.), len(am)))))))
amUpSampled = np.append(am,np.zeros((fs-1,am.size)))
amU = amUpSampled.flatten(0)
#% filtered sequence
st_alpha5 = np.convolve(amU, gt_alpha5)
st_alpha1 = np.convolve(amU, gt_alpha1)
#% taking only the first 10000 samples
st_alpha5 = st_alpha5[0:10000:1]
st_alpha1 = st_alpha1[0:10000:1]
#st_alpha5_reshape = np.reshape(st_alpha5, (fs*2.), (np.dot(N, fs)/20.)).T
st_alpha5_reshape = np.reshape(st_alpha5, (-1,500)).T
#st_alpha1_reshape = np.reshape(st_alpha1, (fs*2.), (np.dot(N, fs)/20.)).T
st_alpha1_reshape = np.reshape(st_alpha1, (-1,500)).T
plt.close('all')
plt.figure(1)
plt.plot(np.array(np.hstack((np.arange(.1, (1.99)+(1./fs), 1./fs)))), np.real(st_alpha5_reshape).T, 'b')
plt.title('eye diagram with alpha=0.5')
plt.xlabel('time')
plt.ylabel('amplitude')
plt.axis(np.array(np.hstack((0., 2., -1.5, 1.5))))
plt.grid('on')
plt.figure(2)
plt.plot(np.array(np.hstack((np.arange(.1, (1.99)+(1./fs), 1./fs)))), np.real(st_alpha1_reshape).T, 'b')
plt.title('eye diagram with alpha=1')
plt.xlabel('time')
plt.ylabel('amplitude')
plt.axis(np.array(np.hstack((0., 2., -1.5, 1.5))))
plt.grid('on')
plt.show()

Please let me know where the issue is and whats the fix is in Python Code only?

Improve this question

1 Answer 1

Reset to default
3

A few things. Despite the fact that you said "I have doubly verified the values of each variables bot in Python and Octave - both of them are same also." -- this is simply not the case.

First, there are times when you need to subtract 1 from your indices when you're porting from MATLAB to numpy, but your code doesn't have any of those.

So everywhere you have something like:

sincOp[int(sincDenZero[0])-1] = 1.

Change it to

sincOp[int(sincDenZero[0])] = 1

Briefly, the reason for this is because the output from np.where is already 0-indexed, so when you subtract 1, you're overcompensating.

Next, you use np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs)))) all over the place, so lets just create a variable and build that once:

fsrange = np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs))))

But that can just be:

fsrange = np.arange(-fs, fs+(1./fs), 1./fs)

Similarly, this huge line:

cosNum = np.cos(np.dot(np.dot(alpha, np.pi), np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs))))))

Can just be:

cosNum = np.cos(alpha * np.pi * fsrange)

And this line:

amUpSampled = np.append(am,np.zeros((fs-1,am.size)))

Should just be (so you don't modify am, and properly specify args to zeros):

amUpSampled = np.vstack([ am, np.zeros([(fs-1.), am.size]) ])

You specified the wrong flatten order here:

amU = amUpSampled.flatten(0)

It should be flattened using FORTRAN order (what MATLAB uses):

amU = amUpSampled.flatten('F')

Same thing when you reshape, you need to specify FORTRAN order:

st_alpha5_reshape = np.reshape(st_alpha5, [(fs*2.), (N * fs / 20.)], 'F').T
st_alpha1_reshape = np.reshape(st_alpha1, [(fs*2.), (N * fs / 20.)], 'F').T

So your corrected python code should look like:

import numpy as np
import matplotlib.pyplot as plt


j = np.complex(0,1)
N = 10**3
#% number of symbols
am = 2.*(np.random.rand(1., N) > 0.5)-1.+np.dot(j, 2.*(np.random.rand(1., N) > 0.5)-1.)
#% generating random binary sequence
fs = 10.
fsrange = np.arange(-fs, fs+(1./fs), 1./fs)

#% sampling frequency in Hz
#% defining the sinc filter
sincNum = np.sin(np.dot(np.pi, fsrange))
#% numerator of the sinc function
sincDen = np.dot(np.pi, fsrange)
#% denominator of the sinc function
sincDenZero = np.where(np.abs(sincDen) < 10**-10);
sincOp=sincNum/sincDen
sincOp[int(sincDenZero[0])] = 1.


#% raised cosine filter
alpha = 0.5
cosNum = np.cos(alpha * np.pi * fsrange)
cosDen = 1.-np.dot(2.*alpha, fsrange)**2.
cosDenZero = np.nonzero(abs(cosDen)<10**-10);
cosOp = cosNum/cosDen
cosOp[int(cosDenZero[0][0])] = np.pi/4.
cosOp[int(cosDenZero[0][1])] = np.pi/4.
gt_alpha5 = sincOp*cosOp

alpha = 1.
cosNum = np.cos(alpha * np.pi * fsrange)
cosDen = 1.-np.dot(2.*alpha, fsrange)**2.
cosDenZero = np.where(abs(cosDen)<10**-10);
cosOp = cosNum/cosDen
cosOp[int(cosDenZero[0][0])] = np.pi/4.
cosOp[int(cosDenZero[0][1])] = np.pi/4.
gt_alpha1 = sincOp*cosOp


#% upsampling the transmit sequence
amUpSampled = np.vstack([ am, np.zeros([(fs-1.), am.size]) ])
amU = amUpSampled.flatten('F')
#% filtered sequence
st_alpha5 = np.convolve(amU, gt_alpha5)
st_alpha1 = np.convolve(amU, gt_alpha1)
#% taking only the first 10000 samples
st_alpha5 = st_alpha5[0:10000]
st_alpha1 = st_alpha1[0:10000]
st_alpha5_reshape = np.reshape(st_alpha5, [(fs*2.), (N * fs / 20.)], 'F').T
st_alpha1_reshape = np.reshape(st_alpha1, [(fs*2.), (N * fs / 20.)], 'F').T
plt.close('all')
X = np.arange(0,1.99, 1.0/fs)
plt.figure(1)
plt.plot(X, np.real(st_alpha5_reshape).T, 'b')
plt.title('eye diagram with alpha=0.5')
plt.xlabel('time')
plt.ylabel('amplitude')
plt.axis(np.array(np.hstack((0., 2., -1.5, 1.5))))
plt.grid('on')

plt.figure(2)
plt.plot(X, np.real(st_alpha1_reshape).T, 'b')
plt.title('eye diagram with alpha=1')
plt.xlabel('time')
plt.ylabel('amplitude')
plt.axis(np.array(np.hstack((0., 2., -1.5, 1.5))))

plt.grid('on')
plt.show()

Which produces the figures you expected.

Side note:

If you have an array/matrix in MATLAB (let's say it's called varname), you can save it to a .mat file with save varname (in MATLAB).

You can then load that array/matrix in python with:

import scipy.io
mat = scipy.io.loadmat("<path of .mat file>")
varname = mat[varname]

You can do this for your entire MATLAB workspace as well with simply save -- in python mat will still just be a dictionary keyed by variable name, so you'd access the individual workspace variables just as you do above.

You can use this to verify, step-by-step what numpy is producing against what you expect it to produce, and figure out what you're doing wrong.

Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

Thanks a lot for the detailed info and guidance. Indeed I learn't so many of my mistakes. Thanks @jedwards a lot.

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.