Matlab画信号语图（时频图）

matlab自带函数spectrogram使用起来不方便，自己尝试画信号语图（时频图）
时频图即每次将信号截取一定长度做傅里叶变换，将该过程重复多次，然后将所得结果组合到一起。

%% 2FSK
clc,clear,close;
fb=1000;%符号速率
fc=25000;%载波频率
fcA=fc-fb;%载波频率1
fcB=fc+fb;%载波频率2
fs=100000;%采样速率
M=fs/fb;%采样倍数
num=100;%产生的符号数量

h_n0=ones(1,M)/M^0.5;%矩形成型
N=(M-1)/2;

%余弦成型
% a=0.8;%滾降系数
% L=13;%保留旁瓣个数
% h_n0=rcosdesign(a,L*2,M,'sqrt');%平方根余弦
% N=L*M;%卷积后的偏移（（滤波器阶数-1）/2)

sigLen=num*M+length(h_n0)-1;%信号长度计算

symb=randi(2,1,num)-1;%产生等概率二进制随机数0、1
sig_A=zeros(1,num*M);
sig_B=zeros(1,num*M);

sig_A(1:M:end)=symb;%产生冲激函数
sig_B(1:M:end)=abs(symb-1);%这里将symb中原本是0的地方转换为1

%产生的符号与矩形进行卷积，即矩形成型
sig_bs_A=conv(sig_A,h_n0);
sig_bs_B=conv(sig_B,h_n0);

crrA=cos((2*pi*fcA*(0:sigLen-1)/fs));%生成载波1
crrB=cos((2*pi*fcB*(0:sigLen-1)/fs));%生成载波2

sig_frqBd=sig_bs_A.*crrA+sig_bs_B.*crrB;%载波与信号相乘，即调制

%显示调制后信号的频谱
sig_frqBd_FFT=fft(sig_frqBd);
sig_FFT_ABS=abs(sig_frqBd_FFT).^2;
fIdx=((-fs/2):(fs/length(sig_FFT_ABS)):(fs/2));
figure;
semilogy(fIdx(1:end-1),fftshift(sig_FFT_ABS),'b-');
title('信号频谱');

%% 2PSK
clc,clear,close;
fb=1000;%符号速率
fc=25000;%载波频率
fs=100000;%采样速率
M=fs/fb;%采样倍数
num=51200;%产生的符号数量

% h_n0=ones(1,M)/M^0.5;%矩形成型
% N=(M-1)/2;

%余弦成型
a=0.8;%滾降系数
L=13;%保留旁瓣个数
h_n0=rcosdesign(a,L*2,M,'sqrt');%平方根余弦
N=L*M;%卷积后的偏移（（滤波器阶数-1）/2)

sigLen=num*M+length(h_n0)-1;

%波形成型
symb=randi(2,1,num)*2-3;%双极性信号
sig=zeros(1,num*M);
sig(1:M:end)=symb;
sig_bs=conv(sig,h_n0);

crr=cos((2*pi*fc*(0:sigLen-1)/fs));%生成载波

sig_frqBd=sig_bs.*crr;%载波与信号相乘，即调制

% 显示调制后信号的频谱
sig_frqBd_FFT=fft(sig_frqBd);
sig_FFT_ABS=abs(sig_frqBd_FFT).^2;
fIdx=((-fs/2):(fs/length(sig_FFT_ABS)):(fs/2));
figure(1);
semilogy(fIdx(1:end-1),fftshift(sig_FFT_ABS),'b-');
title('信号频谱');
% 信号时频图
win_len=512;
%窗函数
w=blackman(win_len).';
%所作次数
spe_times=500;
for spe_loop=1:1:spe_times
    %截取信号的长度
    singnal_temp=sig_frqBd(win_len*(spe_loop-1)+1:win_len*(spe_loop));%sig_frqBd为相应的信号
    fft_spe(spe_loop,:)=20*log10(abs(fftshift(fft(singnal_temp.*w))));
end
figure;
mesh(fft_spe.')
grid on;
view(2)

2FSK信号时频图为：