原创浅析傅里叶频谱分析&滤波&小波分析

 2020-2-13 15:11  4263 23 11 分类: 模拟

浅析傅里叶频谱分析&滤波&小波分析

原始信号由20Hz、50Hz和100Hz三种频率的正弦波组成，如上图所示。

在Figure图上画框截取待分析的数据，如下所示：

浅析傅里叶频谱分析&滤波&小波分析

所得结果如下所示，

浅析傅里叶频谱分析&滤波&小波分析

从上图中可以清晰的看出有三种频率，而且从时频图中可以看到三种频率发生在不同时刻。

点击单边傅里叶频谱，放大频谱如下：

浅析傅里叶频谱分析&滤波&小波分析

画框可以选定滤波范围，如带通滤波，选择如下框，可得结果：

浅析傅里叶频谱分析&滤波&小波分析

MATLAB中具体代码如下：

function frequency_analysis
clc close all
Ts = 0.001;
Fs = 1/Ts;
f1 = 20;
f2 = 50;
f3 = 100;
dt = 0.2;
t1 = (0:Ts:dt-Ts) + 0; t2 = (0:Ts:dt-Ts) + dt;
t3 = (0:Ts:dt-Ts) + 2*dt;
y1 = sin(2*pi*f1*t1);
y2 = sin(2*pi*f2*t2);
y3 = sin(2*pi*f3*t3);
t = [t1 t2 t3];
y = [y1 y2 y3];
figure
plot(t,y)
xlim([t(1) t(end)])
ylim([min(y) max(y)])
xlabel('时间t')
ylabel('信号y(t)')
title('原始信号')
%
% 以下为标准化程序
% 上面的Ts和Fs要给定，后面会用到
set(gcf,'WindowButtonDownFcn',@BtndownFcn);
function BtndownFcn(h,evt)
temppt = get(gca,'CurrentPoint');
startpt.x = temppt(1,1);
startpt.y = temppt(1,2);
endpt = startpt;
height = 0.0001;
width = 0.0001;
rectangle('Position',[startpt.x, startpt.y, width, height],'Tag','rangerect'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);
set(h,'WindowButtonUpFcn',@BtnupFcn);
function BtnmoveFcn(h,evt)
temppt = get(gca,'CurrentPoint');
endpt.x = temppt(1,1);
endpt.y = temppt(1,2);
width = abs(endpt.x-startpt.x)+0.00001;
height = abs(endpt.y-startpt.y)+0.00001;
hrect = findobj('Tag','rangerect');
set(hrect,'Position',[min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);
end
function BtnupFcn(h,evt)
set(h,'WindowButtonMotionFcn','');
set(h,'WindowButtonUpFcn','');
hrect = findobj('Tag','rangerect');
delete(hrect);
BtnUp_Spectrum_Analysis(startpt.x, endpt.x)
end
% 频谱分析
function BtnUp_Spectrum_Analysis(starttime, endtime)
hline = findobj(gca,'type','line');
time = get(hline,'xdata');
laser = get(hline,'ydata');
% 得到截取的分析数据
tempx = find(time >= min(starttime, endtime));
startindex = tempx(1);
tempx = find(time <= max(starttime, endtime));
endindex = tempx(end);
yt = laser(startindex:endindex);
yt = ytmean(yt);
t = time(startindex:endindex);
% 傅里叶变换
[Yf, f] = Spectrum_Calc(yt,Fs);
% 小波变换
scale = 1:50;
cw2 = cwt(yt,scale,'morl');
% 作图
figure
subplot(231)
% 截取的频谱分析的数据
plot(t, yt)
xlim([t(1) t(end)])
ylim([min(yt),max(yt)]);
title('频谱分析数据')
xlabel('时间t')
ylabel('截取的数据y(t)')
h1 = subplot(234);
% 单边傅里叶变换分析频谱
plot(f,Yf,'-')
title('单边傅里叶频谱')
xlabel('频率Hz')
ylabel('|Y(f)|')
set(h1,'ButtonDownFcn',@BtnDown_Filter_fcn)
function BtnDown_Filter_fcn(h,evt)
fig_fft = figure;
plot(f,Yf,'-')
title('单边傅里叶频谱')
xlabel('频率Hz')
ylabel('|Y(f)|')
xlim([0 Fs/2])
ylim([0 max(Yf)])
% 构造右键菜单
filter_flag = 1;
hcmenu = uicontextmenu;
uimenu(hcmenu, 'Label', '低通滤波', 'Callback', @hcb1);
uimenu(hcmenu, 'Label', '高通滤波', 'Callback', @hcb2);
uimenu(hcmenu, 'Label', '带通滤波', 'Callback', @hcb3);
uimenu(hcmenu, 'Label', '带阻滤波', 'Callback', @hcb4);
set(gca,'UIContextMenu',hcmenu);
function hcb1(h,evt)
filter_flag = 1;
end
function hcb2(h,evt)
filter_flag = 2;
end
function hcb3(h,evt)
filter_flag = 3;
end
function hcb4(h,evt)
filter_flag = 4;
end
set(fig_fft,'WindowButtonDownFcn',@BtndownFcn);
function BtndownFcn(h,evt)
if
strcmp(get(h,'SelectionType'),'normal')
temppt = get(gca,'CurrentPoint');
startpt.x = temppt(1,1);
startpt.y = temppt(1,2);
endpt = startpt;
height = 0.0001;
width = 0.0001;
rectangle('Position',
[startpt.x, startpt.y, width, height],'Tag','rangerect_fft'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);
set(h,'WindowButtonUpFcn',@BtnupFcn);
end
function BtnmoveFcn(h,evt)
temppt = get(gca,'CurrentPoint');
endpt.x = temppt(1,1);
endpt.y = temppt(1,2);
width = abs(endpt.xstartpt.x)+0.00001;
height = abs(endpt.ystartpt.y)+0.00001;
hrect = findobj(h,'Tag','rangerect_fft');
set(hrect,'Position', [min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);
end
function BtnupFcn(h,evt)
set(h,'WindowButtonMotionFcn','');
set(h,'WindowButtonUpFcn','');
hrect = findobj(h,'Tag','rangerect_fft');
delete(hrect);
Filter_Analysis(startpt.x, endpt.x);
function Filter_Analysis(x1,x2)
lowfre = min(x1,x2);
highfre = max(x1,x2);
if
lowfre <= 0 || lowfre >= Fs/2
lowfre = 0.01;
end
if
highfre <= 0 || highfre >= Fs/2
highfre= Fs/2-0.01;
end
W1 = lowfre/(Fs/2);
W2 = highfre/(Fs/2);
filter_str = '(低通)';
switch filter_flag
case 1
[b,a] = butter(5,min(W1,W2)); % 低通，画的框的左边为截止频率
filter_str = '(低通)';
case 2
[b,a] = butter(5,max(W1,W2),'high'); % 高通，画的框的右边为截止频率
filter_str = '(高通)';
case 3
[b,a] = butter(5,[W1 W2]); % 带通
filter_str = '(带通)';
case 4
[b,a] = butter(5,[W1 W2],'stop'); % 带阻
filter_str = '(带阻)';
end
y = filter(b,a,yt);
figure
subplot(2,1,1)
plot(t,yt)
xlabel('时间t')
ylabel('信号y')
xlim([t(1) t(end)])
ylim([min(yt),max(yt)]);
title('原始信号')
subplot(2,1,2)
plot(t,y)
xlabel('时间t')
ylabel('信号y')
xlim([t(1) t(end)])
ylim([min(y),max(y)]);
title(sprintf('滤波信号%s%.2fHz~%.2fHz',filter_str,lowfre,highfre))
end
end
end
end
subplot(1,3,[2,3]) % 频率轴化为频率
[X,Y] = meshgrid(t,5/(2*pi)./scale*Fs);
mesh(X,Y,abs(cw2))
view(0,90)
title('时频图')
xlabel('时间')
ylabel('频率')
xlim([t(1) t(end)])
set(gca,'ylim',[0,max(max(Y))])
set(gca,'YScale','log')
set(gca,'YTick', [1:9,10:10:90,100:100:900,1000,2000])
function [Yf, f] = Spectrum_Calc(yt,Fs)
L = length(yt);
NFFT = 2^nextpow2(L);
Yf = fft(yt,NFFT)/L;
Yf = 2*abs(Yf(1:NFFT/2+1));
f = Fs/2*linspace(0,1,NFFT/2+1);
end
end
end
end
其中选框的代码可以标准化，可以用来在用户交互中用户选择范围，代码整理如下
set(gcf,'WindowButtonDownFcn',@BtndownFcn);
function BtndownFcn(h,evt)
temppt = get(gca,'CurrentPoint');
startpt.x = temppt(1,1);
startpt.y = temppt(1,2);
endpt = startpt;
height = 0.0001;
width = 0.0001;
rectangle('Position',[startpt.x, startpt.y, width, height],'Tag','rangerect'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);
set(h,'WindowButtonUpFcn',@BtnupFcn);
function BtnmoveFcn(h,evt)
temppt = get(gca,'CurrentPoint');
endpt.x = temppt(1,1);
endpt.y = temppt(1,2);
width = abs(endpt.x-startpt.x)+0.00001;
height = abs(endpt.y-startpt.y)+0.00001;
hrect = findobj('Tag','rangerect');
set(hrect,'Position',[min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);
end
function BtnupFcn(h,evt)
set(h,'WindowButtonMotionFcn','');
set(h,'WindowButtonUpFcn','');
hrect = findobj('Tag','rangerect');
delete(hrect);
% ProsessFcn(startpt,endpt);
% 选完框后对选择的范围处理
end
end