循环平稳特征检测

function y=cs(x,N,fs,M) %%%%% 参量说明 %%%%%参量怎么取？ % x -- 信号 % N -- 循环谱检测采样长度，必须小于等于信号序列长度 % fs -- 采样频率, 检测带宽为-fs/2至fs/2 % M -- 平滑点数, 时间分辨率*频率分辨率=M %%%%%%%%%%%%%%%%%%% %%% 参数 %%% %N =64; %fs = 16; %M = 64; win = 'hamming'; % 平滑窗类型 if nargin<1;x=randint(1,1024);end; d_alpha = fs/N; % 1/时间分辨率=循环频率分辨率 alpha = 0:d_alpha:30; % 循环频率, 分辨率=1/时间分辨率 a_len = length(alpha); % 循环频率取样个数 f_len = floor(N/M-1)+1; % 最大平滑窗个数, 即频率采样个数 f = -(fs/2-d_alpha*floor(M/2)) + d_alpha*M*(0:f_len-1); % 频率采样点位置 % S = zeros(a_len, f_len); % 初始相关功率谱 S = zeros(f_len,a_len); % 初始相关功率谱 i = 1; %%% 信号fft变换 %%% x = fftshift(fft(x(1:N))); x =x'; %%% 遍历循环频率取值 %%% for alfa = alpha interval_f_N = round(alfa/d_alpha); % 循环频率所对应的频谱序列序号 f_N = floor((N-interval_f_N-M)/M)+1; % 平滑窗的个数 %%% 生成平滑窗函数 %%% g = feval(win, M); %生成M个点的窗 window_M = g(:, ones(f_N,1));%把g的第一列拿出来存f_n列构成新矩阵，f_n为窗的个数 Mxf_n %%% 频域序列平滑模板 %%% t = 1:M*f_N; t = reshape(t, M, f_N); %%% 计算X1,X2 %%% X1 = x(t).*window_M; X2 = x(t+interval_f_N).*window_M; %%% 计算谱相关 %%% St = conj(X1).*X2; St = mean(St, 1); % 平滑平均 % S(i, floor((f_len-f_N)/2)+(1:f_N)) = St/N; % 将结果平移至序列中央以便作图 S(floor((f_len-f_N)/2)+(1:f_N), i) = St/N; % 将结果平移至序列中央以便作图 i = i+1; end %%% 遍历循环频率取值结束 %%% %%% 循环功率谱作图 %%% S=abs(S)./max(max(abs(S))); % mesh(f,alpha,S); figure(4); mesh(alpha,f,S); % mesh(S); % axis tight; %title('BPSK-OFDM'); xlabel('循环频率a'); ylabel('频率f'); title(' 循环谱')