# 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的...