循环平稳 谱切片

李慧聪

请教各位大神！！！！拜谢！！！
  我现在在研究循环平稳（二阶），循环自相关，循环谱，也是在大神们的帮助下做出来了。
但是现在卡在  谱切片   现在快崩溃中，，，希望大神给我指点指点

dsp2008

楼主见过大面包没有？就像切面包一样，对三维循环谱图进行切片。

实际操作一般可这样做，若循环谱图为X(a,f)，其中a是循环频率，f是谱频率，则切片Xc(a)=max(X(:,1:end))。

Raspberry

《基于MATLAB7.0统计信息处理》
第一部分 基础理论
第一章 MATLAB7.0简介
第二章 数理统计基本理论
第三章 高阶统计量理论
第二部分 数理统计工具箱的工程应用
第四章 数理统计工具箱函数简介
第五章 概率分布及其统计特征
第六章 线性统计模型
第七章 非线性统计模型
第八章 假设检验
第九章 多元统计分析
第十章 隐马尔可夫模型
第三部分 高阶统计工具箱的工程应用
第十一章 高阶统计工具箱函数简介
第十二章 高阶统计量的估计
第十三章 非线性相位耦合检测
第十四章 线性过程建模及其应用
第十五章 谐波恢复和波达方向估计
第十六章 非线性系统估计
第十七章 时延估计
第十八章 高阶时频分布及其应用
参考文献
笫12章中有一节仁绍了11/2切片谱的问题！！！！

fearless

切片是要自己根据你的数据手动选择循环频率进行切片，所以进行切片之前，你应该是自己计算好大概会在那个位置切片能够表现出最好的特征。

Triste

问题解决了吗？

think2015

你可以按照二楼说的做一下

李慧聪

谢谢，已经开始做这个了，嘿嘿 ，先用着

李慧聪

dsp2008 发表于 2016-8-15 20:26
楼主见过大面包没有？就像切面包一样，对三维循环谱图进行切片。

实际操作一般可这样做，若循环谱图为X( ...

恩呢 ，我这就试试，谢谢

李慧聪

Raspberry 发表于 2016-8-16 08:40
《基于MATLAB7.0统计信息处理》
第一部分 基础理论
第一章 MATLAB7.0简介

哇 谢大神指教

李慧聪

Triste 发表于 2016-8-17 11:19
问题解决了吗？

没呢，这一段因为一些杂事，瞎忙了。这就开始，谢大神关注

李慧聪

think2015 发表于 2016-8-18 08:27
你可以按照二楼说的做一下

我先试试看嘿嘿  谢谢啦

think2015

李慧聪 发表于 2016-9-8 11:31
我先试试看嘿嘿  谢谢啦

ok  祝你顺利  有问题及时讨论

李慧聪

这个是主程序循环谱密度（但是用这个仿真信号出来的图和正确的图不一样，用其他仿真信号的时候有时又是对的）fs=5120; N_num=4096;  t=[0:1/fs:(N_num-1)/fs];  f1=4.5;f2=8.5;fc=100; x=(1+cos(2.*pi.*f1.*t+cos(2.*pi.*f2.*t)).*cos(2.*pi.*fc.*t));  L = length(x);                % 信号长度 Nw = 256;                        % 窗口长度 Nv = fix(2/3*Nw);        % block overlap 块重叠 nfft = 2*Nw;                % FFT length da = 1/L;           % cyclic frequency resolution循环频率分辨率 a1 = 51;            % first cyclic freq. bin to scan (i.e. cyclic freq. a1*da)                             %起始循环频率 a2 = 500;           % last cyclic freq. bin to scan (i.e. cyclic freq. a2*da)                           %最后的循环频率  % Loop over cyclic frequencies遍历循环频率 C = zeros(nfft,a2-a1+1); S = zeros(nfft,a2-a1+1); Q = ~strcmp(which('chi2inv'),'')==1; % check if function 'chi2inv' is available检查函数的chi2inv是否可用%%strcmp比较字符串 for k = a1:a2;     if Q == 1         Coh = SCoh_W(x,x,k/L,nfft,Nv,Nw,'sym',.01);  %循环谱相关  p=0.1;;p[0,1]     else         Coh = SCoh_W(x,x,k/L,nfft,Nv,Nw,'sym');     end         C(:,k-a1+1) = Coh.C;         S(:,k-a1+1) = Coh.Syx;         waitbar((k-a1+1)/(a2-a1+1))     %？？ end  % Plot results Fs = 5120;                        % sampling frequency in Hz alpha = Fs*(a1:a2)*da; f = Fs*Coh.f(1:nfft/2);  figure imagesc(alpha,f,abs(S(1:nfft/2,:))), colormap(jet),colorbar,axis xy,title('Cyclic Spectral Density'),%循环谱密度 xlabel('cyclic frequency \alpha [Hz]'),ylabel('spectral frequency f [Hz]')  figure imagesc(alpha,f,abs(C(1:nfft/2,:))), colormap(jet),colorbar,axis xy,title('Cyclic Spectral Coherence'),%循环谱相关 xlabel('cyclic frequency \alpha [Hz]'),ylabel('spectral frequency f [Hz]')  figure %ABS=abs(Rx_ALPHA); max_R=max(max(abs(S(1:nfft/2,:)))); contour(abs(alpha),f,abs(S(1:nfft/2,:)),linspace(0.1*max_R,max_R)) ; title('循环谱相关函数{R_{x}}^{\alpha}');xlabel('循环频率 α/Hz');ylabel('频率 τ/s') xlim([0 250]);ylim([0 150]);  figure

李慧聪

子函数   function  Spec = CPS_W(y,x,alpha,nfft,Noverlap,Window,opt,P)   if length(Window) == 1     Window = hanning(Window); end Window = Window(:); n = length(x);          % Number of data points nwind = length(Window); % length of window  % check inputs if (alpha>1||alpha<0),error('alpha must be in [0,1] !!'),end if nwind <= Noverlap,error('Window length must be > Noverlap');end if nfft < nwind,error('Window length must be <= nfft');end if nargin > 7 && (P>=1 || P<=0),error('P must be in ]0,1[ !!'),end  y = y(:); x = x(:);                 K = fix((n-Noverlap)/(nwind-Noverlap));        % Number of windows  % compute CPS index = 1:nwind; f = (0:nfft-1)/nfft; t = (0:n-1)'; CPS = 0; if strcmp(opt,'sym') == 1     y = y.*exp(-1i*pi*alpha*t);     x = x.*exp(1i*pi*alpha*t); else     x = x.*exp(2i*pi*alpha*t); end  for i=1:K     xw = Window.*x(index);     yw = Window.*y(index);     Yw1 = fft(yw,nfft);                % Yw(f+a/2) or Yw(f)     Xw2 = fft(xw,nfft);                % Xw(f-a/2) or Xw(f-a)     CPS = Yw1.*conj(Xw2) + CPS;     index = index + (nwind - Noverlap); end  % normalize KMU = K*norm(Window)^2;        % Normalizing scale factor ==> asymptotically unbiased CPS = CPS/KMU;     % variance reduction factor Window = Window(:)/norm(Window); Delta = nwind - Noverlap; R2w = xcorr(Window); k = nwind+Delta:Delta:min(2*nwind-1,nwind+Delta*(K-1)); if length(k) >1     Var_Reduc = R2w(nwind)^2/K + 2/K*(1-(1:length(k))/K)*(R2w(k).^2); else     Var_Reduc = R2w(nwind)^2/K; end  % confiance interval if nargin > 7     v = 2/Var_Reduc;     alpha = 1 - P;     if alpha == 0                % Sa ~ Chi2         temp = 1./chi2inv([1-alpha/2 alpha/2],round(v));         CI = v*CPS*temp;     else                                    % Sa ~ Normal         Sy = CPS_W(y,y,0,nfft,Noverlap,Window,opt);         Sx = CPS_W(x,x,0,nfft,Noverlap,Window,opt);         Var_Sa = Sy.S.*Sx.S/v;         temp = sqrt(2)*erfinv(2*P-1);         CI = CPS*[1 1] + temp*sqrt(Var_Sa(:))*[1 -1];     end end  % set up output parameters if nargout == 0     figure,newplot;     plot(f,10*log10(abs(CPS))),grid on     xlabel('[Hz]'),title('Spectral Correlation Density (dB)') else     Spec.S = CPS;     Spec.f = f;     Spec.K = K;     Spec.Var_Reduc = Var_Reduc;     if nargin > 7          Spec.CI = CI;     end end

李慧聪

子程序  function Coh = SCoh_W(y,x,alpha,nfft,Noverlap,Window,opt,P)     if length(Window) == 1     Window = hanning(Window); end Window = Window(:); n = length(x); nwind = length(Window);  % check inputs if (alpha>1)||(alpha<0),error('alpha must be in [0,1] !!'),end if nwind <= Noverlap,error('Window length must be > Noverlap');end if nfft < nwind,error('Window length must be <= nfft');end if (nargin>7) && (P>=1 || P<=0),error('P must be in ]0,1[ !!'),end  y = y(:); x = x(:);  t = (0:n-1)'; if strcmp(opt,'sym') == 1     y = y.*exp(-1i*pi*alpha*t);     x = x.*exp(1i*pi*alpha*t); else     x = x.*exp(2i*pi*alpha*t); end  Syx = CPS_W(y,x,0,nfft,Noverlap,Window,opt);     Sy = CPS_W(y,y,0,nfft,Noverlap,Window,opt);      Sx = CPS_W(x,x,0,nfft,Noverlap,Window,opt);      Coh.K = Sx.K;    Coh.f = Sx.f;     Coh.Syx = Syx.S; Coh.Sy = Sy.S; Coh.Sx = Sx.S; Coh.C = Syx.S./sqrt(Sy.S.*Sx.S);  % variance reduction factor Window = Window(:)/norm(Window); Delta = nwind - Noverlap; R2w = xcorr(Window); k = nwind+Delta:Delta:min(2*nwind-1,nwind+Delta*(Coh.K-1)); if length(k) >1     Coh.Var_Reduc = R2w(nwind)^2/Coh.K + 2/Coh.K*(1-(1:length(k))/Coh.K)*(R2w(k).^2); else     Coh.Var_Reduc = R2w(nwind)^2/Coh.K; end  % threshold on square magnitude at P significance level under H0 if nargin > 7     Coh.thres = chi2inv(1-P,2)*Coh.Var_Reduc/2; end  % set up output parameters if nargout == 0     figure,newplot;     subplot(211),plot(Coh.f(1:nfft/2+1),abs(Coh.C(1:nfft/2+1)).^2), grid on     if nargin > 7,         hold on,plot(Coh.f(1:nfft/2+1),Coh.thres,':r'),         title(['Spectral Coherence (squared magnitude) and threshold at ',num2str(100*P),'% level of significance'])     else         title('Spectral Coherence (squared magnitude)')     end     subplot(212),plot(Coh.f(1:nfft/2+1),angle(Coh.C(1:nfft/2+1))), grid on     xlabel('[Hz]'),title('Phase') end