EEPW论坛

fni=input('频域带通滤波－输入数据文件名:','s');
fid=fopen(fni,'r');  
sf=fscanf(fid,'%f',1);%sample frequncy
fmin=fscanf(fid,'%f',1);%最小截止频率（低通0）
fmax=fscanf(fid,'%f',1);%最大截止频率
sx=fscanf(fid,'%s',1);%读入横向坐标的标注（时间/s)
sy=fscanf(fid,'%s',1);%读入纵向坐标的标注（加速度)
fno=facsnf(fid,'%s',1);%输出数据文件名
x=fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
n=length(x);%取信号数据长度
t=(0:1/sf:(n-1)/sf)';%建立离散时间列向量
nfft=2^nextpow2(n);%取大于2并最接近n的2的幂次方为fft长度
ni=round(fmin*nfft/sf+1);%四舍五入取整求最小截止频率对应数组元素下标
na=round(fmax*nfft/sf+1);%四舍五入取整求最大截止频率对应数组元素下标
y=fft(x,nfft);%进行fft变换，结果存于y
a=zeros(1,nfft);
a(ni:na)=y(ni:na);%将y的正频带通内的元素赋值给a
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);%将y的负频带通内的元素赋值给a
y=ifft(a,nfft);%进行fft逆变换，结果存于y
y=(real(y(1:n)))';%取逆变换的实部n个元素为滤波结果列向量
subplot(2,1,1);%绘制滤波前的时程曲线图形
plot(t,x);
xlabel(sx);ylabel(sy);grid on;
subplot(2,1,2);绘制滤波后的时程曲线
plot(t,y);
xlabel(sx);ylabel(sy);grid on;
fid=fopen(fno,'w');
for k=1:n
fprintf(fid,'%f\n',t(k),y(k));
end
status=fclose(fid);