clear;
clc;
close all;
fni=csvread('F:\SRTP\fan data\2016-11-03-不平衡-加配重-2\2016-05-09274-50-1.csv');
fid=fopen(fni,'r');
sf=fscanf(fid,'%f',1); %采样频率
fmin=fscanf(fid,'%f',1); % 最小截止频率
fmax=fscanf(fid,'%d',1); % 最大截止频率
c=fscanf(fid,'%d',1); % 单位变换系数
it=fscanf(fid,'%s',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);% 计算输入数据的长度
nt=length(x);% 建立时间向量
t=0:1/sf:(n-1)/sf;% 取大于n且与其最接近的2的整数次方为FFT长度
nf=2^nextpow2(nt);% FFT 变换
y=fft(x,nfft);% 计算频率间隔
df=sf/nfft;% 计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);% 计算原频率间隔
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=zeros(1,nfft);
a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1);
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);% 积分结果
y=real(y(1:n))*c;
figure;
subplot(211)
plot(t,x);
xlabel(sx);
ylabel(sy1);
grid on;% 积分后曲线
subplot(212)
plot(t,y);
xlabel(sx);
ylabel(sy2);
grid on;% 打开文件输出积分后的数据
fid=fopen(fno,'w');
for k=1:n
fprint(fid,'%f%f\n',t(k),y(k));
end
status=fclose(fid);
%作者:王苏阳
%链接:https://www.zhihu.com/question/37401929/answer/76674187
%来源:知乎
%著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。
评论0