cyclic2_am.m

基于matlab

clear all;clc;clear;
tic %%%计算程序运行时间
% %%% ----------------
% %%% AM 调制
% %%% ----------------
% Fs =10230;
Fs = 5120;
Fc = 200; % Carrier frequency in Hz  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = [1:.1*Fs]'/Fs; % Sampling times for .1 second
x = sin(20*pi*t)'; % Representation of the signal %%%%%%%%%%%%%%%
y = ammod(x,Fc,Fs); % Modulate x to produce y.
figure;
subplot(2,1,1); plot(t,x);title('Representation of the signal'); % Plot x on top.
subplot(2,1,2); plot(t,y);title('Modulate x to produce y'); % Plot y below.

% %-----------------
% %调频
% %-----------------
% clf;
% fm=10;fc=50;fs=1000;N=1024;k=0:N-1;t=k/fs;
% x=sin(2.0*pi*fm*t);
% xf=abs(fft(x,N));
% y=modulate(x,fc,fs,'fm');
% subplot(211);plot(t,y);title('调频信号');
% yf=abs(fft(y,N));
% subplot(212);plot(yf);

% %-----------------
% % 调相
% %-----------------
% clf;
% fm=10;fc=100;fs=1000;N=1024;k=0:N-1;t=k/fs;
% x=sin(2.0*pi*fm*t);
% xf=abs(fft(x,N));
% y=modulate(x,fc,fs,'pm');
% subplot(211);plot(t,y);title('调相信号');
% yf=abs(fft(y,N));
% subplot(212);plot(yf);

% %%---------------
% %%only sezaosheng
% %%---------------
% %N=1024; %采样点数
% N=length(y);
% signalPower=sum(abs(y).^2/length(y));  %信号能量
% % 
% SNR=-8;  
% No=10^(SNR/10);  %单位转化
% noisePower=signalPower/No;   %噪声能量
% noi=sqrt(noisePower^2)*randn(1,N);%白噪声
% %noi=randn(1,N);
% % noi=A*randn(1,N);%高斯白噪声
% 
% p=0.5;%色噪声功率谱谱宽
% f0=0.05; %色噪声功率谱中心频率
% a1=-2*p*cos(2*pi*f0);
% a2=p*p;
% for iii=3:N
%     v(1)=-2;v(2)=-1.8;
%     v(iii)=a1*v(iii-1)-a2*v(iii-2)+noi(iii); %二阶AR模型产生高斯色噪声
% end
% %y=v;
% % figure(1);
% % subplot(211);plot(y); title('色噪声信号的时域图形');
% signal=y+v;

% % %%---------------
% % %% only sin inteference
% % %%---------------
% N=length(y);
% f3=100;%干扰信号频率
% Ay=1; %信号幅度
% SNR=-5;
% No=10^(SNR/20); %单位转化
% An=Ay/No; %干扰信号幅度
% 
% noise=An*sin(2*pi*f3*t); %正弦信号干扰
% % noise=sin(2*pi*f3*t);


signalPower=sum(abs(y).^2/length(y));  %信号能量
N=length(y);

SNR=5;  
No=10^(SNR/10);  %单位转化
noisePower=signalPower/No;   %噪声能量
noise=sqrt(noisePower/2)*randn(1,N);%白噪声
signal=y+noise;
% 
% N=512;
% noise=rand(1,N); %非高斯噪声


save phase signal;


%%%%%%%%%%%%%%%%%%%%%%%%%%
%求循环谱%
%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha_len=(-0.5:1/N:0.5-1/N);
M=N/16; %%%%%%%%%%%%%%%%%  
X=fft(signal);
Y=X;
% figure(5);
% plot(abs(X));
X=fftshift(X);
% figure(2);
% plot(abs(X));title('FFT');

for alpha=1:N/2-M/2+1,
    for f1=1:N/2-M/2+2-alpha,
        if f1==1,
            tmp(alpha,f1)=X([f1+N/2+alpha-1-M/2:f1+N/2+alpha-2+M/2])*(X([f1+N/2-alpha+1-M/2:f1+N/2-alpha+M/2]))';
        else tmp(alpha,f1)=tmp(alpha,f1-1)-X(f1+N/2+alpha-2-M/2)*conj(X(f1+N/2-alpha-M/2))+X(f1+N/2+alpha-2+M/2)*conj(X(f1+N/2-alpha+M/2));
        end;
    end;
end;%tmp为alpha>0,f>0的四分之一平面
% for alpha=1:N/2-M/2+1,
%     for f1=1:N/2-M/2+2-alpha,
%         if f1==1,
%             tmp(alpha,f1)=X([f1+N/2+alpha-1-M/2:f1+N/2+alpha-2+M/2])*(X([f1+N/2-alpha+1-M/2:f1+N/2-alpha+M/2]))';
%         else tmp(alpha,f1)=tmp(alpha,f1-1)-X(f1+N/2+alpha-2-M/2)*conj(X(f1+N/2-alpha-M/2))+X(f1+N/2+alpha-2+M/2)*conj(X(f1+N/2-alpha+M/2));
%         end;
%     end;
% end;%tmp为alpha>0,f>0的四分之一平面

S=[tmp([1:N/2-M/2+1],[N/2-M/2+1:-1:2])/M tmp/M];% f轴扩充
S=[S([N/2-M/2+1:-1:2],:);S]; %alpha轴扩充
S1=zeros(N,N);
S1([M/2:N-M/2],[M/2:N-M/2])=S;
X1=S1(:,N/2+1);
X2=S1(N/2+1,:);
S1=abs(S1);

f_axe=[-0.5:1/N:0.5-1/N]*2*pi;
alpha_axe=(-0.5:1/N:0.5-1/N)*2;
% f_axe=[-2560:5120/N:2560-1/N]*2*pi;
% alpha_axe=(-2560:5120/N:2560-1/N)*2;

figure(3);
subplot(221);contour(f_axe,alpha_axe,S1);xlabel('f'),ylabel('alpha');title('');grid on
subplot(222);mesh(f_axe,alpha_axe,S1);xlabel('f');ylabel('alpha');zlabel('S1');
subplot(223);mesh(f_axe,alpha_axe,S1);xlabel('f');ylabel('alpha');zlabel('S1');view(90,0);
subplot(224);mesh(f_axe,alpha_axe,S1);xlabel('f');ylabel('alpha');zlabel('S1');view(0,0);

figure(4)
subplot(211);plot(alpha_axe,abs(X1));xlabel('alpha');ylabel('abs(X1)');title('');
subplot(212);plot(f_axe,abs(X2));xlabel('f');ylabel('abs(X2)');

toc
% [temp,J]=max(abs(X1));
% X1(J)=0;
% [temp,I]=max(abs(X1));
% xxx=(N/2-I)*2/N
% J
% xxx=(J-I)*1/J