clc;clear all;close all
fs=1000; %%采样频率
f1=50;
f2=100;
Fc=6; %%初始中心频率
w0=2*pi*Fc; %%初始中心频率一般fc=2hz，一把w0>=5,做地震时是w0=2*pi;
dt=1.0/fs;
t=0:dt:1;
Nt=length(t);
nscale=500; %%尺度个数或频率个数
sig = sin(2*pi*f1*t)+sin(2*pi*f2*t);
w_sig=fft(sig); %%第一步，对信号做傅立叶变换
figure(1)
plot(t,sig)
sw_morlet=zeros(Nt,nscale);
smax = Nt*dt; %最大尺度对应最小中心频率
smin = w0*dt/pi; %最小尺度对应最大中心频率 应是要分析的最大频率的2倍
xx=log2(smin); %处理成二进小波尺度 可以减少大量的冗余度
yy=log2(smax);
ss = linspace(xx,yy,nscale);
scale=2.^ss;
Freq=(Fc)./scale; %%此处计算频率轴,也就是横坐标 fc=2，w0=2*pi*fc
% c=2*Fc*nscale;
% scale=c./(1:nscale);
% FreqBins = linspace(smin,smax,nscale);
% scale = Fc./ FreqBins;
% freq = FreqBins * (1/dt);
% scale1=linspace(smin,smax,nscale);
%计算morlet频率小波
for w=0:1:(Nt-1)
for s=1:1:nscale
temp_f=2*pi*w/(Nt*dt); %%计算频率，(0:1:Nt-1)/(Nt*dt)
sw_morlet(w+1,s)=scale(s).^(1/2)*exp(-(scale(s)*temp_f-w0).^2/2);
%sw_morlet(w+1,s)=exp(-(2*scale(s)*pi*w/(Nt*dt)-w0).^2/2);%%计算归一化
end
end
%计算小波系数在频率域
for m=1:1:nscale
c=w_sig.*(sw_morlet(:,m)');
% temp(m,:)=c; %从频率域可以看出两个频率分量在那个尺度而且是否被分开个人觉得看要分开的重点频率是否在滤波器即小波尺度函数的范围内
wcoefs(m,:)=ifft(c);
end
% figure
% a=ifft(sw_morlet(:,20));
% plot(images(a))
[X,Y]=meshgrid(Freq,t);
h=pcolor(X,Y,abs(wcoefs'));
c=jet(256);
colormap(c);
set(h,'LineStyle','none');
box off
% xlim([40 110]);
title('单道时频分析');
set(gca,'ydir','reverse');
% fmax = 0.5;
% fmin = 0.01;
% fb = 4 ;
% fc = 2;
% FreqBins = linspace(fmin,fmax,nscale);
% Scales = fc./ FreqBins;
% RealFreqBins = FreqBins * (1/dt);
% figure(3);
% MWT=cwt(sig,Scales,'cmor4-2');
% pcolor(t,RealFreqBins,abs(MWT));legend('工具箱');colormap jet;shading interp;
% ylabel('Frequency / Hz');
% xlabel('Time / sec');
