[求助]如何获得与50Hz工频同步的窄脉冲信号

[northbear]

如何获得与50Hz工频同步的窄脉冲信号，请高手指点。

[tmaczn]

global dt df t f N
close all

N=2^14; %采样点数
L=32; %每码元的采样点数
M=N/L %码元数
Rb=2; %码速率是2Mb/s
Ts=1/Rb; %码元间隔
dt=Ts/L; %时域采样间隔
df=1/(N*dt) %频域采样间隔
T=N*dt %截短时间
Bs=N*df/2 %系统带宽
Na=4; %示波器扫描宽度为4个码元
alpha=0.5;

t=[-T/2+dt/2:dt:T/2]; %时域横坐标
f=[-Bs+df/2:df:Bs]; %频域横坐标


%升余弦
hr1=sin(pi*t/Ts)./(pi*t/Ts);
hr2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
hr=hr1.*hr2;
HR=abs(t2f(hr));%取模是为了忽略时延
GT=sqrt(HR);
GR=GT;
EP=zeros(size(f));
EP1=EP;

for loop1=1:10
Eb_N0(loop1)=(loop1-1); %Eb/N0 in dB
eb_n0(loop1)=10^(Eb_N0(loop1)/10);
Eb=1;
n0=Eb/eb_n0(loop1); %信道的噪声谱密度
sita=n0*Bs; %信道中噪声功率

n_err=0; %误码计数

for loop2=1:3

a=sign(randn(1,M));
a1=zeros(L,M);
for loop=1:L,a1(loop,:)=a; end
a1=reshape(a1,1,N);

A1=t2f(a1);

P=A1.*conj(A1)/T; %功率谱
EP=(EP*(loop2-1)+P)/loop2; %累计平均
F=30+10*log10(EP+eps); %加eps以避免除以零


%观察信道的眼图
b=sign(randn(1,M))+1;
imp=zeros(1,N); %产生冲激序列
imp(L/2:L:N)=b/dt;

S=t2f(imp).*HR; %升余弦信号的傅氏变换
s=f2t(t2f(imp).*HR); %升余弦信号的时域波形
s=real(s);

%噪声
n_ch=sqrt(sita)*randn(size(t)); %信道噪声
nr=real(f2t(t2f(n_ch).*GR)); %输出噪声

sr=real(f2t(A1.*HR))+nr; %接收信号
y=sr(L/2:L:N); %取样
aa=sign(y); %判决
n_err=n_err+length(find(aa~=a)) ;

y1=zeros(L,M);
for loop=1:L,y1(loop,:)=aa; end
y1=reshape(y1,1,N);
Y1=t2f(y1);
P1=Y1.*conj(Y1)/T; %功率谱
EP1=(EP1*(loop1-1)+P1)/loop1; %累计平均
F1=30+10*log10(EP1+eps); %加eps以避免除以零
end
Pe(loop1)=n_err/(M*loop2);

end





figure(1)
%输入波形
hold on
plot(t,a1,'g')
axis([-7,+7,-2,+2])
xlabel('t(us)')
ylabel('s(t)(V)')


figure(2)
%输入信号为双极性不归零码的功率谱
plot(f,F,'r'); grid
axis([-7,+7,-50,+50])
xlabel('f(MHZ)')
ylabel('Psi(f)(MHZ)')

figure(3)
%输出波形
hold on
plot(t,y1,'g')
axis([-7,+7,-2,+2])
xlabel('t(us)')
ylabel('s(t)(V)')

figure (4) %接收信号的功率谱
plot(f,F1,'r'); grid
axis([-7,+7,-50,+50])
xlabel('f(MHZ)')
ylabel('Pso(f)(MHZ)')

figure(5) %信道的眼图
hold on
tt=[0:dt:4*L*dt];
for jj=1:4*L:N-4*L
plot(tt,s(jj:jj+4*L));
end

figure(6) %输出信号的眼图
hold on
tt=[0:dt:4*L*dt];
for jj=1:4*L:N-4*L
plot(tt,sr(jj:jj+4*L));
end


figure(7)
%判决波形和输出波形的比较
hold on
plot(t,sr,'g')
plot(t,y1,'r');grid
axis([-7,+7,-2,+2])
xlabel('t(us)')
ylabel('s(t)(V)')

figure(8) %误比特率曲线
semilogy(Eb_N0,Pe,'g');
eb_n0=10.^(Eb_N0/10);
hold on
semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)));
axis([0,9,1e-4,1])
xlabel('Eb/N0')
ylabel('Pe')