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

原始信号由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