|
楼主 |
发表于 2012-6-21 09:09
|
显示全部楼层
回复 2 # wanyeqing2003 的帖子
这是我的带通滤波程序
%带通滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('频域带通滤波-输入数据文件名:','s');
fid=fopen(fni,'r'); %fid用于存储文件句柄值
sf=fscanf(fid,'%f',1); %读出文件中的采样频率的信息,读入的是一个1*1的矩阵
fmin=fscanf(fid,'%f',1); %读出文件中的最小截止频率,读入的是一个1*1的矩阵
fmax=fscanf(fid,'%f',1); %读出文件中的最大截止频率,读入的是一个1*1的矩阵
sx=fscanf(fid,'%s',1); %读出文件中的横向坐标轴的标注,读入的是一个1*1的矩阵
sy=fscanf(fid,'%s',1); %读出文件中的纵向坐标轴的标注,读入的是一个1*1的矩阵
fno=fscanf(fid,'%s',1); %输出数据文件名,读入的是一个1*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);
%四舍五入取整求最小截止频率对应数组元素的下标
ni=round(fmin*nfft/sf+1);
%四舍五入取整求最大截止频率对应数组元素的下标
na=round(fmax*nfft/sf+1);
%进行FFT变换,结果存于y
y=fft(x,nfft);
%建立一个长度为nfft元素全为0的向量
a=zeros(1,nfft);
%将y的正频率带通内的元素赋值给a
a(ni:na)=y(ni:na);
%将y的负频率带通内的元素赋值给a
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
%进行fft逆变换,结果存于y
y=ifft(a,nfft);
%取逆变换的实部n个元素为滤波结果列向量
y=(real(y(1: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 %f\n',t(k),y(k));
end
status=fclose(fid);
|
|