% 本程序是用于信号滤波的实验程序。
% 主要功能是进行信号的保相滤波。
% ZhangXining 2001.5.12
% 对计算滤波信号的幅值问题进行了修正
% ZhangXining 2009.12.25
clear all;
pause on;
FileName=input('请输入滤波的数据文件名:','s');
Fs=input('请输入信号的采样频率(Hz):');
load(FileName);
temp=strtok(FileName,'.');
x=eval(temp);
x=x-mean(x);
xlen=max(size(x));
if mod(xlen,2)==1
x=x(1:xlen-1);
xlen=xlen-1;
end
handle1=figure;
t=(0:xlen-1)/Fs;
t=t';
[a b]=size(x);
if a==1
x=x';
end
plot(t, x);
title('振动信号的波形图');
xlabel('时间/s');
ylabel('幅值/mv');
grid on;
disp('按任意键继续......');
pause
%hidem(handle1);
disp('请用鼠标在谱图上选择滤波的截止频率,');
disp('左键选择,右键退出?按任意键继续......');
disp('1. 低通滤波');
disp('2. 高通滤波');
disp('3. 带通滤波');
filtype=input('请选择滤波的类型:');
y=fft(x)/(xlen/2);
absy=abs(y);
angley=angle(y);
z=abs(y(1:xlen/2));
Amp=max(z);
u=ones(1,xlen/2)*Amp;
winlen=round(log2(xlen/40.));
winlen=pow2(winlen);
v=Amp*cos((-pi:(pi/winlen):0))/2+Amp/2;
w=Amp*cos((0:(pi/winlen):pi))/2+Amp/2;
winlen=winlen+1;
f=(0:(xlen-1)/2)*(Fs/xlen);
Df=f(2);
handle2=figure;
plot(f, z);
title('信号的幅值谱图');
xlabel('频率/Hz');
ylabel('幅值/mv');
hold on;
plot(f,u,':m');
r=1;
Fl=0;
Fh=xlen/2;
while r==1
if filtype==3
[X,Y,B]=ginput(1);
XX(1)=X/Df;
BB(1)=B;
[X,Y,B]=ginput(1);
XX(2)=X/Df;
BB(2)=B;
X=XX;
B=BB(1);
if floor(X(1))<((winlen+1)/2)
X(1)=(winlen+1)/2;
end
if ceil(X(1))>(xlen-(winlen+1)/2)
X(1)=xlen-(winlen+1)/2;
end
if floor(X(2))<((winlen+1)/2)
X(2)=(winlen+1)/2;
end
if ceil(X(2))>(xlen-(winlen+1)/2)
X(2)=xlen-(winlen+1)/2;
end
X
else
[X,Y,B]=ginput(1);
X=X/Df;
if floor(X)<((winlen+1)/2)
X=(winlen+1)/2;
end
if ceil(X)>(xlen-(winlen+1)/2)
X=xlen-(winlen+1)/2;
end
end
if B~=1
r=0;
else
if filtype==1 %低通
Fh=ceil(X);
u=ones(1,xlen/2)*Amp;
for i=1:winlen
u(Fh-((winlen-1)/2+1)+i)=w(i);
end
for i=(Fh+(winlen-1)/2+1):xlen/2
u(i)=0;
end
end
if filtype==2 %高通
Fl=ceil(X);
u=ones(1,xlen/2)*Amp;
for i=1:winlen
u(Fl-((winlen-1)/2+1)+i)=v(i);
end
for i=1:(Fl-((winlen-1)/2+1))
u(i)=0;
end
end
if filtype==3 %带通
temp=min(X)
if temp>Fl
Fl=ceil(temp);
for i=1:winlen
u(Fl-((winlen-1)/2+1)+i)=v(i);
end
for i=1:(Fl-((winlen-1)/2+1))
u(i)=0;
end
end
temp=max(X);
if Fh>temp
Fh=ceil(temp);
for i=1:winlen
u(Fh-((winlen-1)/2+1)+i)=w(i);
end
for i=(Fh+((winlen-1)/2+1)):xlen/2
u(i)=0;
end
end
end
end
clf;
plot(f, z);
title('信号的幅值谱图');
xlabel('频率/Hz');
ylabel('幅值/mv');
hold on;
plot(f, u, 'm:');
end
%hidem(handle2);
disp('计算滤波器频谱及脉冲响应函数......');
u=u/Amp;
s=fliplr(u);
s=[u s];
s=ifft(s);
handle3=figure;
plot(real(s(1:100)));
N=input('输入滤波器单位脉冲响应序列的长度:');
%N=N/2;
disp('进行傅里叶反变换,求单位脉冲响应函数');
clf;
s=fliplr(u);
s=[u s];
s=ifft(s);
disps=[s(xlen-N:xlen) s(1:N)];
plot(real(disps),'b');
%hold on;
%plor(imag(disps));
title('滤波器的单位脉冲响应函数');
xlabel('采样点数');
ylabel('幅值/mv');
disp('按任意键继续......');
pause
%hidem(handle3);
disp('按任意键显示滤波结果......');
handle4=figure;
subplot(2,1,1);
plot(t, x, 'g');
ylabel('原始信号');
xlabel('时间/s');
title('滤波结果比较');
grid on;
ft=fliplr(u);
ft=[u,ft];
ft=ft';
ft=ft*(xlen/2);
ft=absy.*ft;
ft=ft.*exp(j*angley);
fs=ifft(ft);
subplot(2,1,2);
plot(t, real(fs), 'b');
ylabel('滤波信号');
xlabel('时间/s');
grid on;
disp('按任意键保存滤波结果......');
pause
%保存滤波后的信号到文件
FiltFile=input('输入保存的文件名:','s');
temp=strtok(FiltFile,'.');
temp=cellstr(temp);
temp=strcat(temp,'.txt');
q=real(fs);
temp=sprintf('save %s %s %s',char(temp),'q','-ascii -double');
eval(temp);
%显示所有的窗口
figure(handle1);
figure(handle2);
figure(handle3);