Преобразование кода Matlab в Python - выход не соответствует

Question

Я успешно преобразовал исходный код MATLAB в Python, но вывод графика просто не соответствует. Я дважды проверил значения каждого бота переменных в Python и Octave - оба они одинаковы.

Октав Plot Выход:

Python Matplotlib Выход:

Октав Код:

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 код:

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() 

Сообщите мне, где проблема и что исправить в Python Code только?

Answer 1

Несколько вещей. Несмотря на то, что вы сказали: «Я дважды проверил значения каждого бота переменных в Python и Octave - оба они одинаковы». -- Это просто не тот случай.

Во-первых, там : раз, когда вам нужно вычесть 1 из ваших индексов при переносе с MATLAB на numpy, но в вашем коде нет ни одного из них.

Так везде у вас есть что-то вроде:

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

изменить его на

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

Вкратце, причина этого в том, что выход из np.where уже 0 индексированные, поэтому при вычитании 1, вы перекомпенсируете.

Далее вы используете np.array(np.hstack((np.arange(-fs, (fs)+(1./fs), 1./fs)))) повсюду, так что позволяет просто создать переменную и строить, что когда-то:

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

Но это может быть просто:

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

Кроме того, эта огромная линия :

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

Может быть просто:

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

И эта линия:

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

Если просто (так что вы не изменяли am и правильно указать арг к zeros):

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

вы указали неправильный расплющить заказ здесь :

amU = amUpSampled.flatten(0) 

Должно быть сплющенные с помощью FORTRAN порядка (что использует MATLAB):

amU = amUpSampled.flatten('F') 

То же самое, когда вы изменить форму, вам нужно указать FORTRAN заказ:

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 

Таким образом, ваш исправленный код питон должен выглядеть следующим образом:

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() 

Который производит цифры, которые вы ожидали.

Примечание стороны:

Если у вас есть массив/матрицы в MATLAB (скажем, это называется varname), вы можете сохранить его в файл .mat с save varname (в MATLAB).

Вы можете загрузить массив/матрицы в питоне с:

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

Вы можете сделать это для всей рабочей области MATLAB, а с просто save - в питоне mat будет еще раз быть словарем ключа по имя переменной, чтобы вы могли обращаться к отдельным переменным рабочей области так же, как и выше.

Вы можете использовать это, чтобы поэтапно проверить, какое количество денег производит против того, что вы ожидаете от него, и выяснить, что вы делаете неправильно.