こんにちは、Python で Matlab コードを変換しようとしています。私のコードは変調までは正常に動作しますが、ノイズを追加してフィルタリングを行うと、間違った結果が得られます。基本的な問題はフィルターにあると思います。
Python コード
from __future__ import division
import numpy as np
from numpy import *
import cmath
from pylab import *
from matplotlib import *
from scipy.signal import firwin ,lfilter, freqz
a1=array([0,0,1])
a2=arange(0,2,1)
N_bits=1e2
bits1=ceil(len(a2)*rand(N_bits))
bits=a1[array(bits1,dtype=int)]
N_bits=len(bits)
delta1=zeros(((N_bits/4)+1,) , dtype=np.complex)
delta2=zeros(((N_bits/4)+1,) , dtype=np.complex)
k=1
C=zeros(((N_bits/4)+1,) , dtype=np.complex)
D=zeros(((N_bits/4)+1,) , dtype=np.complex)
C[0]=2+2j
D[0]=1+1j
for l in arange(0,N_bits,4):
if (bits[l]==0 & bits[l+1]==0):
delta1[k]=0
elif(bits[l]==0 & bits[l+1]==1):
delta1[k]=np.pi/2
elif (bits[l]==1)&(bits[l+1]==1):
delta1[k]=np.pi
else:
delta1[k]=3*(np.pi/2)
if (bits[l+2]==0)&(bits[l+3]==0):
delta2[k]=0
elif(bits[l+2]==0)&(bits[l+3]==1):
delta2[k]=np.pi/2
elif (bits[l+2]==1)& (bits[l+3]==1):
delta2[k]=np.pi;
else:
delta2[k]=3* np.pi/2
C[k]=C[k-1]*(np.cos(delta1[k])+1j*np.sin(delta1[k]))
D[k]=D[k-1]*(np.cos(delta2[k])+1j*np.sin(delta2[k]))
k=k+1
S=C+D
len_S=len(S)
a=int((N_bits/4)+1)
T=20e-3
B=0.99
Fc=8/T
dt=(1/8)/Fc
sa=T/dt
Convfac=5
alph_a=S.real
alph_b=S.imag
alph_aa1=array(zeros((a*sa)))
alph_bb1=array(zeros((a*sa)))
j=0
for k in range(0,len_S):
alph_aa1[j:j+sa-1]=np.kron(ones(1,float,sa),alph_a[k])
alph_bb1[j:j+sa-1]=np.kron(ones(1,float,sa),alph_b[k])
j=j+sa
h_1=alph_aa1.shape
print(len(alph_aa1))
clip_level=1
t1=linspace(-T/2,T/2,sa)
raised=(cos(pi*(t1/T)))**3
f_index=(raised>clip_level).nonzero()
raised[f_index]=clip_level
raised=raised/norm(raised)
alph_aa=(convolve(alph_aa1,raised))/Convfac
alph_aa=alph_aa[(sa/2)-1:-sa/2]
h_2=alph_aa.shape
alph_bb=(convolve(alph_bb1,raised))/Convfac
alph_bb=alph_bb[(sa/2)-1:-sa/2]
print(len(alph_bb))
t=arange(dt,(a*T)+dt,dt)
tt=t.shape
aaa=alph_aa.shape
###############
# Modulation WORKS FINE up to here
mod=alph_aa*cos(2*pi*Fc*t)-alph_bb*sin(2*pi*Fc*t)
phi=(0/180)*pi
Es=10
SNR=10
variance=Es*10**(-SNR/10)
std_dev=sqrt(variance/2)
noise=(np.random.randn(len(mod)))*std_dev
mod_noisy=mod+noise
###############
#Demodulation Problem start here
##############
demod_1_a=mod_noisy*2*cos(2*pi*Fc*t+phi)
N=10
Fc=40
Fs=1600
# Filter PART consider this part
d=firwin(numtaps=N,cutoff=40,nyq=Fs/2)
Hd=lfilter( d, 1.0, demod_1_a)
y2=(convolve(Hd,raised))/Convfac
y2=y2[(sa/2)-1:-sa/2]
demod_3_a=y2[(sa/2)-1:sa:,]
demod_1_b=-1*mod_noisy*2*sin(2*pi*Fc*t+phi)
Hd=lfilter(d,1.0,demod_1_b)
y2=(convolve(Hd,raised))/Convfac
y2=y2[(sa/2)-1:-sa/2]
demod_3_b=y2[(sa/2)-1:sa:,]
#########3333
#Demod
demod=demod_3_a+(1j)*demod_3_b
print(demod)
Matlab コード:
%%%%%%%%%%%%%%ENCODING%%%%%%%%%%%%%%%%
a1=[0,1];
N_bits=1e2;
bits=a1(ceil(length(a1)*rand(1,N_bits)));
delt1=zeros(1,N_bits/4);
delt2=zeros(1,N_bits/4);
k=2;
C=zeros(1,N_bits/4);
D=zeros(1,N_bits/4);
C(1)=2+2j;
D(1)=1+1j;
for i=1:4:N_bits
if bits(i)==0&&bits(i+1)==0
delt1(k)=0;
elseif bits(i)==0&&bits(i+1)==1
delt1(k)=pi/2;
elseif bits(i)==1&&bits(i+1)==1
delt1(k)=pi;
else
delt1(k)=3*pi/2;
end
if bits(i+2)==0&&bits(i+3)==0
delt2(k)=0;
elseif bits(i+2)==0&&bits(i+3)==1
delt2(k)=pi/2;
elseif bits(i+2)==1&&bits(i+3)==1
delt2(k)=pi;
else%if x(3)==1&&x(4)==0
delt2(k)=3*pi/2;
end
C(k)=C(k-1)*(cos(delt1(k))+1j*sin(delt1(k)));
D(k)=D(k-1)*(cos(delt2(k))+1j*sin(delt2(k)));
k=k+1;
end
S=C+D;
% S=S(1,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=N_bits/4+1;
T=20e-3;
B=0.99;
Fc=8/T;
dt=1/8/Fc;
sa=T/dt;
ConvFac=5
alph_a=real(S);
alph_b=imag(S);
aplh_aa1=zeros(1,a*sa);
alph_bb1=zeros(1,a*sa);
j=1;
for k=1:a
alph_aa1(1,j:j+sa-1)=repmat(alph_a(k),1,sa);
alph_bb1(1,j:j+sa-1)=repmat(alph_b(k),1,sa);
j=j+sa;
end
h_1=size(alph_aa1)
%%%%%%%%% raised cosine %%%%%%%%%
clip_level=1;
t1=linspace(-T/2,T/2,sa);
raised=cos(pi*t1/T).^3
f_index= find(raised>clip_level)
raised(f_index)=clip_level
raised=raised/norm(raised);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%convolving with raised cosine%%%%%
alph_aa=conv(alph_aa1,raised)/ConvFac;
alph_aa=alph_aa(sa/2:end-(sa/2));
h_2=size(alph_aa)
alph_bb=conv(alph_bb1,raised)/ConvFac;
alph_bb=alph_bb(sa/2:end-(sa/2));
%%%%%%%%%%%modulation%%%%%%%%%%%
t=dt:dt:a*T;
length(t)
tt=size(t)
aaa=size(alph_aa)
mod=alph_aa.*cos(2*pi*Fc.*t)-alph_bb.*sin(2*pi*Fc.*t);
%%%%%%%%%%%%%%%%%%%CHANNEL%%%%%%%%%%%%%
phi=0/180*pi;
%mod=awgn(mod,10);
Es=10;
SNR=10;
variance=Es*10^(-SNR/10);
std_dev=sqrt(variance/2);
noise=(randn(1,length(mod)))*std_dev;
mod_noisy=mod+noise;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%demodulation%%%%%%%%%%%%
demod_1_a=mod_noisy*2.*cos(2*pi*Fc*t+phi);
d=fdesign.lowpass('N,Fc',10,40,1600);
Hd = design(d);
y = filter(Hd,demod_1_a);
y2=conv(y,raised)/ConvFac;
y2=y2(sa/2:end-(sa/2));
demod_3_a=y2(sa/2:sa:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%
demod_1_b=-1*mod_noisy*2.*sin(2*pi*Fc*t+phi);
y = filter(Hd,demod_1_b);
y2=conv(y,raised)/ConvFac;
y2=y2(sa/2:end-(sa/2));
demod_3_b=y2(sa/2:sa:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
demod=demod_3_a+1i*demod_3_b
フィルターが適切に設計されていないと思います。復調とフィルタリングの観点からコードを検討してください。