modest.m

有关于相位

clear all;
close all;
clc;
fc=2000; %载波频率
fs=60000; %采样速率
ts=1/fs;                               %采样周期
t0=5.5;                               %信号长度
t=[0:ts:t0];

fd=125;                              %符号速率
ti=1/fd;                               %码元间隔
M=32;                                %码元个数
N=ti/ts;                            %周期比

select=menu('调制方式','2ASK','2FSK','2PSK','4ASK','4FSK','4PSK', '16QAM');
switch select
case 1,                                % 2ASK signal
x=randint(1,M);                     %随机码元
m=sin(2*pi*fc*t);                   %载波

y=zeros(1,M*N);
for i=1:M
    for j=1:N
       y((i-1)*N+j)=x(i)*m(j);
    end
end

case 2,                                 %2FSK signal
x=randint(1,M);  
m1=sin(2*pi*fc*t);
m2=sin(2*pi*2*fc*t);

y=zeros(1,M*N);
for i=1:M
   if x(i)==1;
    for j=1:N;
        y((i-1)*N+j)=m1(j);
    end
   elseif  x(i)==0;
     for j=1:N;
       y((i-1)*N+j)=m2(j);
     end 
  end
end

case 3,                                 %2PSK signal,
x=randint(1,M); 
m1=sin(2*pi*fc*t);
m2=sin(2*pi*fc*t+pi);
y=zeros(1,M*N);
for i=1:M
   if x(i)==1;
    for j=1:N;
        y((i-1)*N+j)=m1(j);
   end
   elseif  x(i)==0;
     for j=1:N;
       y((i-1)*N+j)=m2(j);
     end 
   end
end

case 4,                                 % 4ASK signal
x=randint(1,M, 4);
m=sin(2*pi*fc*t);
y=ones(1,M*N);
for i=1:M
    for j=1:N
       y((i-1)*N+j)=x(i)*m(j);
  end  
end

case 5,                                 % 4FSK signal
x=randint(1,M, 4);
m1=sin(2*pi*fc*t);
m2=sin(2*pi*2*fc*t);
m3=sin(2*pi*3*fc*t);
m4=sin(2*pi*4*fc*t);
y=zeros(1,M*N);
for i=1:M
    if x(i)==0;
    for j=1:N
       y((i-1)*N+j)=m1(j);
    end
elseif x(i)==1;
    for j=1:N
       y((i-1)*N+j)=m2(j);
    end
elseif x(i)==2;
    for j=1:N
       y((i-1)*N+j)=m3(j);
    end
elseif x(i)==3;
    for j=1:N
       y((i-1)*N+j)=m4(j);
    end  
  end  
end

case 6,                                     %4PSK signal
x=randint(1,M, 4);
m1=sin(2*pi*fc*t);
m2=sin(2*pi*fc*t+pi/2);
m3=sin(2*pi*fc*t+pi);
m4=sin(2*pi*fc*t+3*pi/2);
y=zeros(1,M*N);
for i=1:M
    if x(i)==0;
    for j=1:N
       y((i-1)*N+j)=m1(j);
    end
elseif x(i)==1;
    for j=1:N
       y((i-1)*N+j)=m2(j);
    end
elseif x(i)==2;
    for j=1:N
       y((i-1)*N+j)=m3(j);
    end
elseif x(i)==3;
    for j=1:N
       y((i-1)*N+j)=m4(j);
    end  
  end  
end;

case 7,                                     %16QAM signal
x1 = randint(1, M, 4);
x1(find(x1==0)) = -1;
x1(find(x1==2)) = -3;

x2 = randint(1, M, 4);
x2(find(x2==0)) = -1;
x2(find(x2==2)) = -3;

m1=sin(2*pi*fc*t);
m2=sin(2*pi*fc*t+pi/2);
y=zeros(1,M*N);

for i=1:M
    for j=1:N
       y((i-1)*N+j)=x1(i)*m1(j) + x2(i)*m2(j);
    end
end

end
SNR=10;                                    %定义信噪比，单位DB 
sigpow=mean(abs(y).^2);                       %power of input signal
noisefac=10^(-SNR/10);
noise=randn(1,size(y,2));
noise=noise*(sqrt(sigpow*noisefac)/sqrt(mean(noise.^2)));      %产生所需的高斯噪声
ynoise=noise+y;                               %加噪后的信号

% y = ynoise;

hxn = hilbert(y);
% z = sqrt(mod .* mod + imag(hxn).^2);
% %plot(z);
% 
% ma = mean(z);
% an = (z./ma-1).^2;
% %plot(an);
% 
% NUM = length(z);
% y1 = fft(an)/NUM;
% 
% max(abs(y1))
% std(an)
% phase = unwrap(angle(hxn))-2*pi*fc/fs;
rawPhase = unwrap(angle(hxn));
% 
rawFreqV = diff(rawPhase);
meanFreq = mean(rawFreqV);

index = find(abs(rawFreqV) < 1.05 * abs(meanFreq) & abs(rawFreqV) > 0.95 * abs(meanFreq));
freqV = rawFreqV(index);
phase = rawPhase(index);

FreqWndLen = N * 2;
freqWnd = ones(1, FreqWndLen) * meanFreq;
lastPT = 1;
crPhase = 0;

% plot(rawFreqV);
% figure;
% plot(freqV);

for n = 1 : length(freqV) - FreqWndLen
%     if ( abs(freqV(n)) < 1.1 * abs(meanFreq) & abs(freqV(n)) > 0.9 * abs(meanFreq))
%         meanFreq = (meanFreq * FreqWndLen - freqWnd(lastPT) + freqV(n)) / FreqWndLen; 
%         freqV(n) - freqWnd(lastPT)
%         FreqWnd(lastPT) = freqV(n);
%         lastPT = lastPT + 1;
%         lastPT = mod(lastPT, FreqWndLen);
%         if (lastPT == 0)
%             lastPT = FreqWndLen;
%         end
% %     end; 
   
    meanFreq= mean(freqV(n: n + FreqWndLen));
%     meanFreqV(n) = meanFreq;
    crPhase(n+1) = crPhase(n) + meanFreq; 
end
% crPhase(end) = [];
phase(length(crPhase)+1:end) = [];
% 
% % plot(freqV); grid on;
zz = phase - crPhase;
% plot(zz);grid on;

% plot(1:length(freqV), phase, 1:length(freqV), crPhase);


% phase1 = phase(481:end);
% phase2 = phase(1:end-480);
% 
% phase3 = phase2 - phase1;
% 
% zz = phase3;
zz1 = floor(zz/2/pi);
zz2 = zz - zz1 * 2*pi;
zz2(find(zz2>2*pi-0.1)) = 0;

plot(zz2);grid on;

% freqV = diff(phase);
% % freq = mean(freqV);
% 
% index = find(abs(freqV)<pi/4);
% phase1 = phase(index);
% freqV1 = freqV(index);
% freq1 = mean(freqV1);
% % yy1 = y(index);
% Num = length(index);
% % freq1=(phase(end) - phase(1))/Num;
% yy(1) = 0;
% for i=2:Num
%     yy(i) = yy(i-1) + freq1;
% end;
% zz = phase1 - yy;
% zz1 = floor(zz/2/pi);
% zz2 = zz - zz1 * 2*pi;
% 
% % plot(zz);grid on;
% % figure;
% plot(zz2);grid on;
% % figure;
% % plot(1:Num, phase, 1:Num, yy);grid on;
% 
% ma = sum(zz2.^2)/Num