声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1610|回复: 1

[综合讨论] matlab频域微分程序

[复制链接]
发表于 2014-12-28 22:24 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我有一段matlab求频域微分的程序不太懂,求大神指点一二,不胜感激
%频域微分
%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

clc

close all hidden
format long

%%%%%%%%%%%%%%%%%%%%%%

fni=input('频域微分-输入数据文件名:','s');

fid=fopen(fni,'r');

sf = fscanf(fid,'%f',1);     %采样频率

fmin = fscanf(fid,'%f',1);   %最小截止频率
fmax = fscanf(fid,'%f',1);   %最大截止频率

c = fscanf(fid,'%f',1);      %单位变换系数

it = fscanf(fid,'%f',1);     %微分次数

sx = fscanf(fid,'%s',1);     %横向坐标轴的标注   

sy1 = fscanf(fid,'%s',1);    %纵向坐标轴输入单位的标注   

sy2 = fscanf(fid,'%s',1);    %纵向坐标轴输出单位的标注   

fno = fscanf(fid,'%s',1);    %输出数据文件名

x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量

status=fclose(fid);

n=length(x);

%建立时间向量

t=0:1/sf:(n-1)/sf;            

%大于并最接近n的2的幂次方为FFT长度   

nfft=2^nextpow2(n);           
%FFT变换

y=fft(x,nfft);               
%计算频率间隔(Hz/s)

df=sf/nfft;
%计算指定频带对应频率数组的下标

ni=round(fmin/df+1);      

na=round(fmax/df+1);
%计算圆频率间隔(rad/s)   

dw=2*pi*df;
%建立从负到正的圆频率向量

w1=0:dw:2*pi*(0.5*sf-df);

w2=2*pi*(0.5*sf-df):-dw:0;

w=[w1,w2];
%以微分次数为指数,建立圆频率变量向量   

w=w.^it;                     
%进行微分的频域变换

a=y.*w;                        

if it==2                                                

%进行二次微分的相位变换   

   y=-a;                       

else                          

%进行一次微分的相位变换                              
  real(y)=imag(a);            
  imag(y)=-real(a);

end

a=zeros(1,nfft);
%消除指定正频带外的频率成分                       

a(ni:na)=y(ni:na);            
%消除指定负频带外的频率成分                  

a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
%IFFT变换                    

y=ifft(a,nfft);   
%取逆变换的实部n个元素并乘以单位变换系数为微分结果  
y=real(y(1:n))*c;            
%绘制微分前的时程曲线图形                    

subplot(2,1,1);                                                            
plot(t,x);                    
%添加横向坐标轴的标注

xlabel(sx);                  
%添加纵向坐标轴的标注   

ylabel(sy1);                  

grid on;                     
%绘制微分后的时程曲线图形                        
subplot(2,1,2);                                                                 
plot(t,y);  
%添加横向坐标轴的标注

xlabel(sx);                  
%添加纵向坐标轴的标注   

ylabel(sy2);                  

grid on;                     
%打开文件输出微分后的数据                        
fid=fopen(fno,'w');           

for k=1:n

fprintf(fid,'%f %f\n',t(k),y(k));

end

status=fclose(fid);

回复
分享到:

使用道具 举报

发表于 2014-12-29 10:27 | 显示全部楼层
每行都有注释啊
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-28 20:50 , Processed in 0.077547 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表