matlab开发-spectral

  • j4_951220
    了解作者
  • 2.6KB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • VIP专享
    资源类型
  • 0
    下载次数
  • 2022-04-02 10:01
    上传日期
matlab开发-spectral。周期图法频谱
spectral.zip
  • spectral.m
    6.4KB
内容介绍
function varargout = spectral(x,dt,win,Nb,P,Wn,ftype,n) % Spectrum using periodogram method % % USAGE: % q = spectral(x,dt,win,Wn,ftype,n) % [psdf,f] = spectral(x,dt,win,Wn,ftype,n) % [psdf,psdfc,f] = spectral(x,dt,win,Wn,ftype,n) % % DESCRIPTION: % Calculates the spectrum for x % using the periodogram method % If a window other than boxcar is used % than the method is refered to as modified % periodogram method. % The confidence intervals are calculated % using the inverse of chi-square CDF. % Also includes a filtering option using the % butterworth filter to see the effect of the % filter on the spectrum % % INPUT VARIABLES: % x - Time series, [vector] % dt - Sampling Rate, [scalar] % win - Window, one of: % 'hanning', 'hamming', 'boxcar' % Nb - Band Averaging, number of bands to average % P - Probability for confidence intervals % Wn - Cut-Off frequencies, used for filtering % ftype - Type of filter, 'high', 'low' or 'stop' % ncb - Number of coefficients to use in % the Butterworth filter % % OUTPUT VARIABLES: % q - structure with the following fields: % xp - detrended x % f = Frequencies % T - Periods % m - Magnitude % a - Amplitude % s - Power spectrum, Sxx(win), [Power] % psdw - Power Spectral Density, Pxx(win), [Power/rad/sample] % psdf - Power Spectral Density, Pxx(f), [Power/sample-freq] % psdT - Power Spectral Density, Pxx(T), [Power*time-unit] % conf - Upper and Lower Confidence Interval multiplication % factors using chi-squared approach % %Copy-Left, Alejandro Sanchez if nargin<2 dt = 1; end if nargin<3 win = 'boxcar'; end if nargin<4 || isempty(Nb) Nb = 1; end if nargin<5 || isempty(P) P = 0.95; end if (nargin<7) && (nargin>5) ftype = 'low'; end if (nargin<8) && (nargin>5) n = 5; end if nargin==0 spectraldemo return end x = x(:); isf = isfinite(x); if nargout==0 figure t = dt*(0:(length(x)-1)); subplot(4,1,1) plot(t,x,'b') hold on t = t(isf); end x = x(isf); mx = mean(x); x = detrend(x); %also removes mean N = sum(isf); fn = (1/(2*dt)); %Nyquist frequency nfft = max(256, 2^(nextpow2(N))); f = (1:nfft/2)/(nfft/2)*fn; %frequencies, Note: starts from 1 f = f(:); %Frequencies %Makes a window try w = eval([win,'(',num2str(N),')']); catch error('Invalid Window') end %Butterworth Filter if nargin>5 [b1,b2] = butter(n,Wn*2*dt,ftype); %Wn/nyquist = Wn*2*dt x = filtfilt(b1,b2,x); xf = x; end if nargout==0 plot(t,x+mx,'r') axis tight set(gca,'xlim',[t(1),t(end)]) end q.xp = x; %save before applying window x = w.*x; %apply window % Evaluate the window normalization constant % A 1/nfft factor has been omitted since it will cancel below R = w'*w; % compute sum of squares of window %Compute Fourier Coefficients %Pad with zeros to length nfft, makes fft much faster %and increases frequency resolution, df=1/N/dt; c = fft(x,nfft); %The magnitude is multiplied by one of the following constants %to compensate for the lost energy when the window is applied %the values are taken from Emery & Thomson pag. 446-448 %In the case of the boxcar window the constant is 1 and is %therefore ommitted if strcmpi(win,'hanning') K = 1/2*sqrt(8/3); elseif strcmpi(win,'hamming') K = 1/2*sqrt(5/2); elseif strcmpi(win,'boxcar') K = 1; else error('Invalid Window') end %Magnitude m = K*abs(c(2:nfft/2+1)); %skip first coefficient which is the mean %Do band averaging now if Nb>1 f = bandaverage(f,Nb); m = bandaverage(m,Nb); end T = 1./f; %Periods %Amplitude a = m.*2./R; %Sxx(win), Angular-Frequency Power Spectrum, [Power] s = 2*(m.^2)./R; %Correction for one-sided spectrum s(2:end-1) = 2*s(2:end-1); %Pxx(win), Amgular-Frequency Spectrum %Power Spectrum Density (PSD), [Power/rad/sample] psdw = s/2/pi; %Pxx(f), Frequency Spectrum %Power Spectrum Density (PSD), [Power/sample-freq] psdf = s*dt; %Pxx(T), Period Spectrum %Power Spectrum Density (PSD), [Power*time-unit] psdT = 2*pi*psdw./T.^2; %Finally Calculate phase ph = unwrap(angle(c(1:fix(N/2)))); %Confidence Intervals using chi-squared approach v = 2*Nb; alpha2 = (1-P)/2; conf = v./chi2inv([alpha2,1-alpha2],v); %Collect variables switch nargout case 0 %do nothing case 1 varargout{1}.f = f; varargout{1}.T = T; varargout{1}.m = m; varargout{1}.a = a; varargout{1}.s = s; varargout{1}.psdw = psdw; varargout{1}.psdf = psdf; varargout{1}.psdT = psdT; varargout{1}.conf = conf; case 2 varargout{1} = psdf; varargout{2} = f; case 3 varargout{1} = psdf; varargout{2} = conf; varargout{3} = f; otherwise error('Invalid number of output variables') end %switch %Plotting if nargout==0 subplot(4,1,2) semilogx(T,a) set(gca,'xlim',[T(end),T(1)]) set(gca,'xticklabel',get(gca,'xtick')) xlabel('Period') ylabel('Amplitude') subplot(4,1,3) loglog(T,psdT) hold on loglog(T,psdT*conf,'--r') axis tight set(gca,'xticklabel',get(gca,'xtick')) xlabel('Period') ylabel('PSD(T)') %Filter or window subplot(4,1,4) if exist('Wn','var') H = freqz(b1,b2,f*pi*2*dt); semilogx(T,abs(H),'b'); axis tight set(gca,'xticklabel',get(gca,'xtick')) xlabel('Period') ylabel('Filter Response') else plot(w) ylabel('window') set(gca,'xticklabel',get(gca,'xtick')) axis tight % plot(ph) % xlabel('Period') end end return %---------------------------------------------------- function yy = bandaverage(y,Nb) %Does the band averaging Lb = floor(length(y)/Nb); yy = reshape(y(1:Nb*Lb),Nb,Lb); yy = mean(yy,1); y = y(:); return %------------------------------------------- function spectraldemo dt = 1/100; t=0:dt:200; sr = 2*randn(size(t)); %white noise s1 = 4*sin(2*pi*t/1); s2 = 2*sin(2*pi*t/4); s3 = 1*sin(t*2*pi/10); s4 = 3*sin(t*2*pi/20); s5 = -0.05*t; s0 = 5; x = s0 + s1 + s2 + s3 + s4 + s5 + sr; x(1000:6000) = NaN; Wn = 1; ftype = 'low'; win = 'hanning'; Nb = []; P = 0.95; n = 5; spectral(x,dt,win,Nb,P,Wn,ftype,n) return
评论
    相关推荐
    • Matlab合集
      冈萨雷斯数字图像处理MATLAB版.中文版+数字图像处理第二版中文版(冈萨雷斯)+MATLAB-R2014a完全自学一本通+MATLAB R2016a完全自学一本通 素材文件+[模式识别与智能计算:MATLAB技术实现(第2版)].杨淑莹.扫描版
    • MATLAB教程
      MATLAB教程MATLAB教程MATLAB教程MATLAB教程MATLAB教程MATLAB教程
    • MATLAB
      MATLAB 该项目是在matlab上完成的,涉及创建和移动宇宙飞船和机器人。 太空飞船和机器人是使用Matlab中的简单几何形状创建的,并通过连续变换矩阵进行移动。 这个项目教我如何使用变换矩阵(旋转,平移等)的概念...
    • MATLAB基础
      一本学习matlab的一本好书
    • MATLAB编译器
      基于MATLAB 2018b版本介绍MATLAB编译器。介绍如何利用编译器将MATLAB代码编译为独立应用程序或组件,并在没有安装MATLAB的计算机上进行部署。
    • matlabruntime
      通过奇点容器运行您的matlab项目 可以在没有MATLAB的容器中运行matlab代码。 为此,我们首先需要通过Matlab编译器在本地构建相应的Matlab代码的独立应用程序,然后使用具有Matlab运行时( 的容器)运行该应用程序 )...
    • matlab实现
      matlab实现 matlab实现matlab实现matlab实现matlab实现
    • matlab 教程
      matlab 信号处理资料,里面包含信号处理pdf文档,一些杂乱的程序和命令等
    • matlab教程
      matlab教程,PPT格式,包含 matlab基本知识、matlab入门、matlab作图、线性规划、无约束优化、非线性规划、统计工具箱、差值、微分方程等多项知识点,并且每个知识点独立成为PPT,内还含有matlab信号处理详解等文档...
    • matlab简介
      1.MATLAB 开发环境 1.1 MATLAB 的视窗环境 进入MATLAB之后,会看到一个视窗MATLAB Command Window称为指令视窗,它是你键入指令的地方同时 MATLAB也将计算结果显示在此。 1.2 简易计算 我们先从MATLAB的...