循环普的 估计fam

FAM的matlab的实现 可以对循环普进行估计 来计算循环普 和其他方法做比较

function [alpha0alpha_profile]=fam（FsdfX）

%Sx=fam（FsdfdalphaNpLPNX）;

% computes spectral auto correlation （SCD） of Signal X

% X input
% Fs sampling frequency
% df Desired frequency resolution
% dalpha Desired cyclic frequency resolution

% SX output
% alpha Cyclic Frequency
% fo Spectrum Frequency

%-------------ver description-------------%
% input SSB modulation also added.
%-----------------------------------------%



% clear all;clc; close all;
% load lp_512Hz_fs8192Hz.dat;
% Fs=8192;
% df=256;
% dalpha=64;
Np=pow2（nextpow2（Fs/df））;
L=Np/4; % This should be chosen to be less than or equal to Np/4
N=length（X）;
N=pow2（nextpow2（N））;
P=pow2（nextpow2（N/L））;
%N=P*L;
% 
% % Insert input signal here%
% %%----------------------------------------------------%%
% fm=512;
% fc=2048;
% A_dc=.5;
% 
% %  sig=sin（2*pi*（fm/Fs）.*（1:N））‘;
%   c=cos（2*pi*（fc/Fs）.*（1:N））‘;
% %  sig_phaseshift=cos（2*pi*（fm/Fs）.*（1:N））‘;
% %  c_phaseshift=sin（2*pi*（fc/Fs）.*（1:N））‘;
% 
%  x=randn（N1）;
%  sig=filter（lp_512Hz_fs8192Hz1x）;
% %AM_sig=（A_dc+sig）.*c; %%%%%%%%%%ordinary AM Modulation
% AM_sig=sig.*c;%%%%%%%%%DSB-SC
% %AM_sig=sig.*c+sig_phaseshift.*c_phaseshift; %%%%%%%%%SSB
% %%----------------------------------------------------%%
AM_sig=X;
if length（AM_sig）    AM_sig（N）=0;
elseif length（AM_sig）>N
    AM_sig=AM_sig（1:N）;
end
NN=（P-1）*L+Np;
xx=AM_sig;
xx（NN）=0;
xx=xx（:）;
X=zeros（NpP）;

for k=0:P-1
    X（:k+1）=xx（k*L+1:k*L+Np）;
end

% windowing ----------------------------------%
w=hamming（Np）;
%XW=diag（w）*X;
XW=X;

%-----------------First FFT--------------------%
XF1=fft（XW）;
XF1=fftshift（XF1）;
XF1=[XF1（:P/2+1:P） XF1（:1:P/2）];

%-----------------Down conversion--------------%
E=zeros（NpP）;
i=sqrt（-1）;
for k=-Np/2:Np/2-1
    for m=0:P-1
        E（k+Np/2+1m+1）=exp（-i*2*pi*k*m*L/Np）;
    end
end
XD=XF1.*E;
XD=conj（XD‘）;
clear （‘XF1‘‘E‘‘X‘）;

%----------------Multiplication----------------------------%

XM = zeros（PNp^2）;
for k = 1:Np
    for c = 1:Np
        XM（:（k-1）*Np+c） = （XD（:k）.*conj（XD（:c）））;
    end
end
%---------------------------------------------------------%

%-------------------Second FFT----------------------------%
XF2 = fft（XM）;
XF2 = fftshift（XF2）;
XF2 = [XF2（:Np^2/2+1:Np^2） XF2（:1:Np^2/2）];
XF2 = XF2 （P/4:3*P/4:）;
MM = abs（XF2）;

alpha0 = -Fs:Fs/N:Fs;
f0 = -Fs/2:Fs/Np:Fs/2;
Sx = zeros（Np+1 2*N+1）;

for k1 = 1:P/2+1
    for k2 = 1:Np^2
        if rem（k2Np） == 0
            c = Np/2 - 1;
        else
            c = rem（k2Np） - Np/2 - 1;
        end
        k = ceil（k2/Np）-Np/2-1;
        p = k1-P/4-1;
        alpha = （k-c）/Np+（p-1）/N;
        f = （k+c）/2/Np;
        if alpha<-1 | alpha>1
            k2 = k2+1;
        elseif f<-.5 | f>.5
            k2 = k2+1;
        else
            kk = 1+Np*（f + .5）;
            ll = 1+N*（alpha + 1）;
            Sx（round（kk） round（ll）） = MM（k1k2）;
        e

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4134  2009-06-03 22:08  fam.m