ask.m

主要用来进行MATLAB中AWGN的识别算法!

% ask.m
% Simulation program to build BPSK modulation singal
%************************ Preparation part *****************************
NUM = 1024;            % The number of FFT point
sr=1000000.0;          % Symbol rate
fc = 1000000;          % The frequency of the carrier
fs = 32000000;         % The sample rate
ml=1;                  % Number of modulation levels
br=sr.*ml;             % Bit rate (=symbol rate in this case)
nd=NUM*fc/fs;         % The number of symbols
n =fc/sr;              % The number of the carrier in a symbol
ebno=14;               % Eb/N0
IPOINT=32;             % Number of oversamples
Ts = 1/sr;             % The period of the sourse signal

t=0:1/fs:1/fc*nd*n-1/fs;    % The sample point in one carrier period
fcos = cos(2*pi*fc*t);      % cosin carrier

%*********************** Filter initialization **************************

irfn=21;                                 % Number of filter taps
alfs=0.33;                                % Rolloff factor
[xh]=hrollfcoef(irfn,IPOINT,sr,alfs,1);  % Transmitter filter coefficients

%*********************** Data generation **********************************
data=rand(1,nd)>0.5;                     % rand: built in function
%*********************** Bpsk modulation **********************************

[data2]=oversample(data,nd,IPOINT);
data3=conv(data2,xh);                                %(conv: build in function)
data4=data3(irfn*IPOINT/2-2:NUM+irfn*IPOINT/2-3);    % date4 is the baseband singal of BPSK 
modulation = data4.*fcos;                            % The modulation is the BPSK singal

%*********************** Attenuation Calculation **************************

spow=sum(modulation.*modulation)/nd;
attn=0.5*spow*sr/br*10.^(-ebno/10);
attn=sqrt(attn);

%*********************** Add White Gaussian Noise(AWGN) *******************

inoise=randn(1,length(modulation)).*attn;
modulation=modulation+inoise;

%fprintf('%f,0,\n',modulation);

%*********************** compute the ampitude  **********************************

modulation1=modulation+0.001;

%fprintf('%f,0,\n',modulation1);
hxn=hilbert(modulation);
y=imag(hxn);
x=real(hxn);
z=sqrt(modulation.*modulation+y.*y);                 %z is ampitute

%****************compute the rmax(零中心归一化瞬时幅度之谱密度的最大值）************************

a=sum(z);
ma=a/NUM;

for i=1:length(z) 
  an(i)=z(i)./ma;
  an(i)=an(i)-1;
  an(i)=an(i)*an(i);
end
  %fprintf('%f\n',an);
  y1=fft(an,NUM);
  yxk=abs(y1(1:NUM/2));
  
 for i=1:NUM/2
    for j=1:NUM/2-i
        t=yxk(j);yxk(j)=yxk(j+1);yxk(j+1)=t;
    end
end

max=10*yxk(NUM/2)/NUM;
fprintf('\nrmax(零中心归一化瞬时幅度之谱密度的最大值)=%f',max);
if(max>0.5)
    fprintf('\nthe signal is not 2FSK,continue to recognize');

%*********************** modulation * modulation **********************************

modulation=modulation.*modulation;
modulation=modulation-0.5;                          %0.5 is the  DC   
hxn=hilbert(modulation);
y=imag(hxn);                               %y is the hilbert change
x=real(hxn);                               %x is the signal ask  

%*********************** compute the phase **********************************
zph=atan(y./x);
z_ph1=y./x;
for i=1:NUM
    if (y(i)>0)&(x(i)>0)
        z_ph(i)=atan(z_ph1(i));
    elseif (y(i)>0)&(x(i)<0)
        z_ph(i)=pi+atan(z_ph1(i));
    elseif (y(i)<0)&(x(i)<0)
        z_ph(i)=pi+atan(z_ph1(i));
    elseif (y(i)<0)&(x(i)>0)
        z_ph(i)=2*pi+atan(z_ph1(i));
    elseif (y(i)>0)&(x(i)==0)
        z_ph(i)=0.5*pi;
    elseif (y(i)<0)&(x(i)==0)
        z_ph(i)=1.5*pi;
    end                                        %  z_ph is the phase
end

%*********************** compute the ddp(零中心瞬时相位的标准差）*******************************

fai=sum(z_ph);
fai0=fai/NUM;                                  %fai0 is the avg
z_ph=z_ph-fai0;  
ddp=sqrt(sum(z_ph.*z_ph)/NUM-sum(z_ph)/NUM*sum(z_ph)/NUM);
fprintf('\nddp(零中心瞬时相位的标准差)=%f',ddp);
if(ddp<1.5)
    fprintf('\n the signal is 2ASK');
end
end
%end the program