马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
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);
|