精华内容
下载资源
问答
  • 短时频谱分析,频谱图的伪彩色映射及显示.频谱图的类型有宽带和窄带两种,伪彩色显示的映射 可按需求设置.文中介绍了基于FFT频谱分析与频谱图显示的原理,Matlab的相关功能函数,并 给出了一个程序实例及其实验结果....
  • %%%%%%傅里叶变换/逆变换/短时傅里叶变换%%%%%%%[y,Fs]=wavread('C:\Users\HSF\Desktop\AgilentTechnologies.wav'); %读出信号,采样率和采样位数。[y,Fs]=wavread('C:\Users\HSF\Desktop\sound\beiyou.wav');%读出...

    %%%%%%傅里叶变换/逆变换/短时傅里叶变换%%%%%%

    %[y,Fs]=wavread('C:\Users\HSF\Desktop\Agilent

    Technologies.wav'); %读出信号,采样率和采样位数。

    [y,Fs]=wavread('C:\Users\HSF\Desktop\sound\beiyou.wav');

    %读出信号,采样率和采样位数。

    sound(y,Fs);

    y=y(:,1);

    %我这里假设你的声音是双声道,我只取单声道作分析,如果你想分析另外一个声道,请改成y=y(:,2)

    figure(1)

    sigLength=length(y);

    t=(0:length(y)-1)/Fs; %计算时间轴

    subplot(3,1,1);plot(t,y);xlabel('Time(s)');%在第一个窗口画波形

    grid on;

    n=0:sigLength-1;

    Y= fft(y,length(y));

    %在sigLength这个有限区间内做快速傅立叶变换

    mag=abs(Y);

    f=Fs*n/sigLength;

    subplot(3,1,2);plot(f,abs(Y));

    xlabel('Frequency(Hz)');

    %在第一个窗口画率谱,

    grid on;

    subplot(3,1,3);

    xifft=ifft(Y);

    sound(xifft,Fs);

    plot(t,xifft);xlabel('Time(s)')

    grid on;

    figure(2)

    Nw=20;

    %窗函数长

    window length

    L=Nw/2;

    %窗函数每次移动的样点数

    Ts=round((sigLength-Nw)/L)+1;

    %计算把数据y共分成多少段

    nfft=512;

    % FFT的长度

    TF=zeros(Ts,nfft);

    %将存放三维图谱,先清零

    for i=1:Ts

    xw=y((i-1)*L+1:i*L+L);   %取一段数据

    temp=fft(xw,nfft);

    %

    FFT变换

    %temp=fftshift(temp);

    %频谱以0频为中心

    TF(i,:)=temp;

    %把谱图存放在TF中

    代表TF矩阵中的第i行

    end

    fn=(1:nfft)*Fs/nfft;

    tn=(1:Ts)*Nw/2/Fs;

    [T,F]=meshgrid(tn,fn);

    mesh(T,F,abs(TF.'));  %三维绘图

    axis tight;

    title('三维时频图');

    ylabel('Frequency(Hz)');

    xlabel('Time(s)');

    grid on;

    仿真结果图如下:​

    a4c26d1e5885305701be709a3d33442f.png

    a4c26d1e5885305701be709a3d33442f.png

    展开全文
  • Matlab中,做短时傅里叶变换需要使用函数spectrogram,而在Matlab2019中,引入了一个新的函数stft,下面我们就来看下这两个函数都如何使用。 短时傅里叶变换的基本原理就是将数据分段加窗,做fft,在分段时会有...

    05b59153bd25b9ce1f6ec1dbc7ee37a1.png

    在Matlab中,做短时傅里叶变换需要使用函数spectrogram,而在Matlab2019中,引入了一个新的函数stft,下面我们就来看下这两个函数都如何使用。

    短时傅里叶变换的基本原理就是将数据分段加窗,做fft,在分段时会有overlap,因此一个向量的短时傅里叶变换结果是一个矩阵。了解了这点,下面的函数及参数就更加容易理解了。

    spectrogram

    参数列表

    先来看spectrogram函数,在更早期的版本中,这个函数的名字是specgram,几种常用的用法如下:

    spectrogram(x)
    s = spectrogram(x)
    s = spectrogram(x, window)
    s = spectrogram(x, window, noverlap)
    s = spectrogram(x, window, noverlap, nfft)
    s = spectrogram(x, window, noverlap, nfft, fs)
    [s, f, t] = spectrogram(x, window, noverlap, nfft, fs)
    [s, f, t] = spectrogram(x, window, noverlap, f, fs)
    [s, f, t, p] = spectrogram(x, window, noverlap, f, fs)
    

    其中,

    • x表示输入信号;
    • window表示窗函数,如果window的值是一个整数,那么被分段的x的每一段的长度都等于window,并采用默认的Hamming窗;如果window是一个向量,那么被分段后每一段的长度都等于length(window),且输入的向量即为所要加的窗函数;
    • overlap表示两段之间的重合点数,overlap的值必须要小于窗长,如果没有指定overlap,默认是窗长的一半,即50%的overlap;
    • nfft表示fft的点数,fft的点数跟窗长可以是不同的,当没有指定该参数时,Matlab会取max(256, 2^(ceil(log2(length(window))))),即当窗长小于256时,fft的点数是256;当窗长大于256时,fft的点数取大于窗长的最小的2的整数次幂;
    • fs表示采样率,用来归一化显示使用;
    • f表示显示的频谱范围,f是一个向量,长度跟s的行数相同;
      • 当x是实信号且nfft为偶数时,s的行数为(nfft/2+1)
      • 当x是实信号且nfft为奇数时,s的行数为(nfft+1)/2
      • 当x是复信号时,s的行数为nfft
      • 当在输入的参数列表中指定f后,函数会在f指定的频率处计算频谱图,返回的f跟输入的f是相同的;
    • t表示显示的时间范围,是一个向量,长度跟s的列数相同;
    • p表示功率谱密度,对于实信号,p是各段PSD的单边周期估计;对于复信号,当指定F频率向量时,P为双边PSD;如何计算PSD

    Examples

    首先,生成信号如下,4个点频信号拼接起来:

    clc;clear all;close all;
    fs = 10e6;
    n = 10000;
    f1 = 10e3; f2 = 50e3; f3 = 80e3; f4 = 100e3;
    t = (0:n-1)'/fs;
    sig1 = cos(2*pi*f1*t);
    sig2 = cos(2*pi*f2*t);
    sig3 = cos(2*pi*f3*t);
    sig4 = cos(2*pi*f4*t);
    
    sig = [sig1; sig2; sig3; sig4];
    

    信号的时域波形如下:

    09afcb7363cd684f060cce033b7cd0d9.png

    直接调用spectrogram(sig),可得如下结果,图中默认横轴是频率,纵轴是时间

    df0d0a17bfe581c9d48be7557576a472.png

    为了绘图更灵活,我们不直接用spectrogram绘图,而且求出s后,再对s单独绘图,这次我们指定window的大小为256

    s = spectrogram(sig, 256);
    t = linspace(0, 4*n/fs, size(s,1));
    f = linspace(0, fs/2, size(s,2));
    figure;
    imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
    colorbar;
    

    2cd0e8dfe447a844f570bed32bf9f31d.png

    noverlap默认是50%,现在我们把它设为window的长度减1,即每次的步进为1

    s = spectrogram(sig, 256, 255);
    t = linspace(0, 4*n/fs, size(s,1));
    f = linspace(0, fs/2, size(s,2));
    figure;
    imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
    colorbar;
    

    58a2e0e60b9feb78fc892b72cfa8fcfc.png

    再加上nfftfs参数,我们指定fft点数就是窗长

    s = spectrogram(sig, 256, 128, 256, fs);
    

    这个的图形跟之前一样,不再画了

    如果在返回值中,增加ft,这样我们下面就不用再重新定义ft

    [s, f, t] = spectrogram(sig, 256, 128, 256, fs);
    figure;
    imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
    colorbar;
    

    从上面的图中我们可以看出,我们的4个信号的频率都比较小,而画出来的图显示的频谱范围比较大,导致下面很大一部分信息我们其实都不需要。这时,我们就可以通过指定f的区间来计算频谱。为了显示效果更好,我们把其他参数也调一下

    window = 2048;
    noverlap = window/2;
    f_len = window/2 + 1;
    f = linspace(0, 150e3, f_len);
    [s, f, t] = spectrogram(sig, window, noverlap , f, fs);
    figure;
    imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
    colorbar;
    

    bd5ddb60f7d2129d19121a8169271e1f.png

    最后再把功率谱密度的返回值加上

    [s, f, t, p] = spectrogram(sig, window, nfft, f, fs);
    figure;
    imagesc(t, f, p);xlabel('Samples'); ylabel('Freqency');
    colorbar;
    

    344f02a3badd551321964a33345a1088.png

    stft

    这个函数在Matlab的解释并不是很多,example也只写了两个,但用法比较简单:

    window = 2048;
    noverlap = window/2;
    nfft = window;
    [s, f, t, p] = spectrogram(sig, window, noverlap, nfft, fs);
    figure;
    imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
    title('使用spectrogram画出的短时傅里叶变换图形');
    colorbar;
    
    ss = stft(sig,fs,'Window',hamming(window),'OverlapLength',window/2,'FFTLength',nfft);
    figure;
    imagesc(t, f, 20*log10((abs(ss(1024:end,:)))));xlabel('Samples'); ylabel('Freqency');
    title('使用stft画出的短时傅里叶变换图形');
    colorbar;
    

    a73e2df027a0197b68f68b483cd2d785.png

    欢迎关注微信公众号:Quant_Times

    http://weixin.qq.com/r/Mi7t9T3EpWDarXmW93sg (二维码自动识别)

    展开全文
  • 实现FFT点数可调、步进可调的短时傅里叶变换;实现加窗功能,窗类型可选。 代码如下 Fs=1000; %采样率1KHz Ts=1/Fs; %点时间间隔0.001s L=2048; %2048点数 t=(0:L-1)Ts; %采的时间为20480.001=2.048s %t1=(0:L/4-13)...

    实现FFT点数可调、步进可调的短时傅里叶变换;实现加窗功能,窗类型可选。
    代码如下
    Fs=1000; %采样率1KHz
    Ts=1/Fs; %点时间间隔0.001s
    L=2048; %2048点数
    t=(0:L-1)Ts; %采的时间为20480.001=2.048s
    %t1=(0:L/4-13)Ts; %500个点
    %t2=(0:L/4+35)Ts; %548个点
    t1=(0:L/4-13)Ts; %512-13=499;500个点
    t2=(0:L/4+35)Ts; %512+35=547;548个点
    %y=sin(50
    2
    pi
    t)+1; %信号周期:T1=50Hz,T2=120Hz;幅度范围为-2~2;
    %y=sin(50
    2pit)+sin(1202pit)+2; %输出无符号十进制,范围0~4
    %x = [2
    cos(2 * pi * 10 * t1+pi/2),2cos(2 * pi * 20 * t1+pi/2),2cos(2 * pi * 30 * t1+pi/2),2cos(2 * pi * 40 * t1+pi/2)]+2;%输出无符号十进制,范围0~4;2048个点
    x = [2
    sin(2 * pi * 10 * t1),2sin(2 * pi * 20 * t1),2sin(2 * pi * 30 * t1),2sin(2 * pi * 40 * t2)]+2;
    y = x
    1023/8192;
    sinwave=round((y2^12));%输出有符号十进制数据;12位数据,最大2047;幅度-2~2
    %%%%%%%%%%循环FFT
    PointNum=1024;
    Steps=256;
    Times=5;
    Waveform_Original_Record=zeros(1,Times
    PointNum);
    Waveform_Windowing_Record=zeros(1,TimesPointNum);
    fft_phase_arranged=zeros(1,Times
    PointNum);
    fftr_amplitude_arranged=zeros(1,TimesPointNum);
    f=zeros(1,Times
    PointNum);
    %加窗处理
    Hanning = hann(PointNum,‘symmetric’);
    Hs_Trans=[Hanning.’ Hanning.’];
    %r1=round(Hs_Trans.sinwave);
    for iteration = 0:1:Times-1
    Waveform_Original_Record(1,iteration
    PointNum+1:iterationPointNum+PointNum)=sinwave(1,iterationSteps+1:iterationSteps+PointNum);
    Waveform_Windowing=(Hanning.’).sinwave(1,iterationSteps+1:iteration
    Steps+PointNum);
    Waveform_Windowing_Record(1,iterationPointNum+1:iterationPointNum+PointNum)=Waveform_Windowing;
    fft_result=fft(Waveform_Windowing); %fft运算
    fft_phase=angle(fft_result);
    fft_phase_arranged(1,iterationPointNum+1:iterationPointNum+PointNum)=[fft_phase(PointNum/2+1:PointNum) fft_phase(1:PointNum/2)];
    fftr_amplitude=abs(fft_result/PointNum); %幅度归一,双边累加
    fftr_amplitude_arranged(1,iterationPointNum+1:iterationPointNum+PointNum)=[fftr_amplitude(PointNum/2+1:PointNum) fftr_amplitude(1:PointNum/2)]; %取正半边
    %f_pos=(1:L/4)Fs/(L/2);
    %f_inv=-1
    fliplr(f_pos);
    %f=[f_inv f_pos];
    f_pos=(1:PointNum/2);
    f_inv=-1fliplr(f_pos);
    f(1,iteration
    PointNum+1:iterationPointNum+PointNum)=[f_inv f_pos]+PointNum/2+iterationPointNum;
    end
    %%%%%%%%%单次FFT
    % a=7;
    % r1=(Hs.’).sinwave(1,a128+1:a128+1024);
    % fftr=fft(r1); %fft运算
    % fftr_normalized=fftr/1024;
    % fft_divide=angle(fftr);
    % fft_phase=[fft_divide(L/4+1:L/2) fft_divide(1:L/4)];
    % %fftr_mod=2
    abs(fftr2/L); %幅度归一,双边累加
    % fftr_mod=2
    abs(fftr2/L); %幅度归一,双边累加
    % fft_pos=[fftr_mod(L/4+1:L/2) fftr_mod(1:L/4)]; %取正半边
    % %f_pos=(1:L/4)Fs/(L/2);
    % %f_inv=-1
    fliplr(f_pos);
    % %f=[f_inv f_pos];
    % f_pos=(1:L/4);
    % f_inv=-1
    fliplr(f_pos);
    % f=[f_inv f_pos];
    %%%%%%%%%%%%%计算公式
    %f=(1:L/4)Fs/(L/2); %频率信息解析,此时步进为Fs/(L/2),点数为L/4。
    %r=ceil(y
    (2^12-1)/4); %数据位宽:12bit
    %%%%%%%%%%%%%%%%%%%%%%%%
    fid = fopen(‘wave_sin.coe’,‘w’);%输出coe格式文件
    fprintf(fid,‘MEMORY_INITIALIZATION_RADIX=10;\n’);
    fprintf(fid,‘MEMORY_INITIALIZATION_VECTOR=\n’);
    for i = 1:1:2^11
    if sinwave(i)>=0
    fprintf(fid,’%d’,sinwave(i));
    else
    fprintf(fid,’%d’,0);
    end
    if i==2^11
    fprintf(fid,’;’);
    else
    fprintf(fid,’,’);
    end
    if mod(i,15)==0
    fprintf(fid,’\n’);
    end
    end
    fclose(fid);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % subplot(3,2,1);
    % plot(1:L/2,sinwave(1,1:L/2));
    % xlabel(‘Point’);
    % ylabel(‘Amplitude’)
    % title(“原始波形”,‘FontSize’,16);
    % subplot(3,2,2);
    % plot(1:L/2,r1(1,1:L/2));
    % xlabel(‘Point’);
    % ylabel(‘Amplitude’)
    % title(“加窗后波形”,‘FontSize’,16);
    % subplot(3,2,3);
    % plot(f(1:1024),fft_pos(1:1024));
    % xlabel(‘Point’);
    % ylabel(‘Amplitude’)
    % title(“一次FFT幅度数据”,‘FontSize’,16);
    % subplot(3,2,4);
    % plot(f(1:1024),fft_phase(1:1024));
    % xlabel(‘Point’);
    % ylabel(‘Phase’)
    % title(“一次FFT相位数据”,‘FontSize’,16);
    % subplot(2,2,1);
    % plot(f,Waveform_Original_Record);
    % xlabel(‘Point’);
    % ylabel(‘Amplitude’)
    % title(“原始波形”,‘FontSize’,16);
    % subplot(2,2,2);
    % plot(f,Waveform_Windowing_Record);
    % xlabel(‘Point’);
    % ylabel(‘Amplitude’)
    % title(“加窗后波形”,‘FontSize’,16);
    subplot(2,2,1);
    plot(f(1:PointNum),fftr_amplitude_arranged(1:PointNum));
    xlabel(‘Point’);
    ylabel(‘Amplitude’)
    title(“首次FFT幅度数据”,‘FontSize’,16);
    subplot(2,2,2);
    %plot(f(12PointNum+1:12PointNum+PointNum),fft_phase_arranged(12PointNum+1:12PointNum+PointNum));
    plot(f(1:PointNum),fft_phase_arranged(1:PointNum));
    xlabel(‘Point’);
    ylabel(‘Phase’)
    title(“首次FFT相位数据”,‘FontSize’,16);
    subplot(2,2,3);
    plot(f,fftr_amplitude_arranged);
    xlabel(‘Point’);
    ylabel(‘Amplitude’)
    title(“一轮STFT幅度数据”,‘FontSize’,16);
    subplot(2,2,4);
    %plot(f(12PointNum+1:12PointNum+PointNum),fft_phase_arranged(12PointNum+1:12PointNum+PointNum));
    plot(f,fft_phase_arranged);
    xlabel(‘Point’);
    ylabel(‘Phase’)
    title(“一轮STFT相位数据”,‘FontSize’,16);
    % subplot(2,2,4);
    % plot(1:L/2,fftr);

    效果图如下
    在这里插入图片描述

    展开全文
  • 首先根据他这个代码和我之前手上已经拥有...短时傅里叶变换,其实还是傅里叶变换,只不过把一段长信号按信号长度(nsc)、重叠点数(nov)重新采样。 % 结合之前两个版本的stft,实现自己的周期图法,力求通俗...

    参考文章:http://www.cnblogs.com/adgk07/p/9314892.html

    首先根据他这个代码和我之前手上已经拥有的那个代码,编写了一个适合自己的代码。

    首先模仿他的代码,测试成功。

    思路:

    短时傅里叶变换,其实还是傅里叶变换,只不过把一段长信号按信号长度(nsc)、重叠点数(nov)重新采样。

    % 结合之前两个版本的stft,实现自己的周期图法,力求通俗易懂,代码分明。
    % 该代码写的时候是按照输入信号为实数的思路写的,在每个片段fft时进行前一半行的转置存储。后续代码思路已修改。 clear; clc; close all
    % %---------------------Chirp 信号生成(start)--------------------- fa = [ 50 800 ]; % 扫描频率上下限 fs = 6400; % 采样频率 % 分辨率相关设定参数 te = 1; % [s] 扫频时间长度 fre = 8; % [s] 频率分辨率 tre = 0.002; % [Hz] 时间分辨率 t = 0:1/fs:te; % [s] 时间序列 sc = chirp(t,fa(1),te,fa(2)); % 信号生成 %待传参(实参) signal = sc; nsc = floor(fs/fre); nff = 2^nextpow2(nsc);%每个窗口进行fft的长度 nov = floor(nsc*0.9); % %---------------------Chirp 信号生成(end)--------------------- %% 使用fft实现周期图法 % Input % signal - Signal vector 输入信号(行向量) % nsc - Number SeCtion 每个小分段的信号长度 % nov - Number OverLap 重叠点数 % fs - Frequency of Sample 采样率 % Output % S - A matrix that each colum is a FFT for time of nsc % 如果nfft为偶数,则S的行数为(nfft/2+1),如果nfft为奇数,则行数为(nfft+1)/2 % 每一列是一个小窗口的FFT结果,因为matlab的FFT结果是对称的,只需要一半 % W - A vector labeling frequency % T - A vector labeling time % Signal Preprocessing h = hamming(nsc, 'periodic'); % Hamming weight function 海明窗加权,窗口大小为nsc L = length(signal); % Length of Signal 信号长度 nstep = nsc-nov; % Number of STep per colum 每个窗口的步进 ncol = fix( (L-nsc)/nstep ) + 1; % Number of CoLum 信号被分成了多少个片段 nfft = 2^nextpow2(nsc); % Number of Fast Fourier Transformation 每个窗口FFT长度 nrow = nfft/2+1; %Symmetric results of FFT STFT_X = zeros(nrow,ncol); %STFT_X is a matrix,初始化最终结果 index=1;%当前片段第一个信号位置在原始信号中的索引 for i=1:ncol %提取当前片段信号值,并用海明窗进行加权(均为长度为nsc的行向量) temp_S=signal(index:index+nsc-1).*h'; %对长度为nsc的片段进行nfft点FFT变换 temp_X=fft(temp_S,nfft); %从长度为nfft点(行向量)的fft变换中取一半(转换为列向量),存储到矩阵的列向量 STFT_X(:,i)=temp_X(1:nrow)'; %将索引后移 index=index + nstep; end % Turn into DFT in dB STFT1 = abs(STFT_X/nfft); STFT1(2:end-1,:) = 2*STFT1(2:end-1,:); % Add the value of the other half STFT3 = 20*log10(STFT1); % Turn sound pressure into dB level % Axis Generating fax =(0:(nfft/2)) * fs/nfft; % Frequency axis setting tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs; % Time axis generating [ffax,ttax] = meshgrid(tax,fax); % Generating grid of figure % Output W = ffax; T = ttax; S = STFT3; %% 画频谱图 subplot(1,3,1) % 绘制自编函数结果 my_pcolor(W,T,S) caxis([-130.86,-13.667]); title('自编函数'); xlabel('时间 second'); ylabel('频率 Hz'); subplot(1,3,2) % 绘制 Matlab 函数结果 s = spectrogram(signal,hamming(nsc),nov,nff,fs,'yaxis'); % Turn into DFT in dB s1 = abs(s/nff); s2 = s1(1:nff/2+1,:); % Symmetric results of FFT s2(2:end-1,:) = 2*s2(2:end-1,:); % Add the value of the other half s3 = 20*log10(s2); % Turn sound pressure into dB level my_pcolor(W,T,s3 ) caxis([-130.86,-13.667]); title('Matlab 自带函数'); xlabel('时间 second'); ylabel('频率 Hz'); subplot(1,3,3) % 两者误差 my_pcolor(W,T,20*log10(abs(10.^(s3/20)-10.^(S/20)))) caxis([-180,-13.667]); title('error'); ylabel('频率 Hz'); xlabel('时间 second'); suptitle('Spectrogram 实现与比较');

    内部调用了my_pcolor.m

    function [  ] = my_pcolor( f , t , s )
    %MY_PCOLOR 绘图函数
    % Input            
    %       f               - 频率轴矩阵
    %       t               - 时间轴矩阵
    %       s               - 幅度轴矩阵
        gca = pcolor(f,t,s);                      % 绘制伪彩色图
            set(gca, 'LineStyle','none');         % 取消网格,避免一片黑
        handl = colorbar;                         % 彩图坐标尺
            set(handl, 'FontName', 'Times New Roman', 'FontSize', 8)
            ylabel(handl, 'Magnitude, dB')        % 坐标尺名称
    end
    View Code

     

    function [  ] = my_pcolor( f , t , s )
    %MY_PCOLOR 绘图函数
    % Input            
    %       f               - 频率轴矩阵
    %       t               - 时间轴矩阵
    %       s               - 幅度轴矩阵
        gca = pcolor(f,t,s);                      % 绘制伪彩色图
            set(gca, 'LineStyle','none');         % 取消网格,避免一片黑
        handl = colorbar;                         % 彩图坐标尺
            set(handl, 'FontName', 'Times New Roman', 'FontSize', 8)
            ylabel(handl, 'Magnitude, dB')        % 坐标尺名称
    end

     

    接下来在主体不变的情况下,修改输入信号和函数返回结果,使其和自己之前的代码效果相同。

     

    其实只要搞清楚了理论,实现上很简单。

    一、首先分析Matlab周期图法API。

    [spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);

    如果函数没有返回结果的调用,MATLAB直接帮我们绘图了。

    1)当输入信号为复数时的代码:

    clear
    clc
    close all
    Fs = 1000;            % Sampling frequency
    T = 1/Fs;             % Sampling period
    L = 1000;             % Length of signal
    t = (0:L-1)*T;        % Time vector
    S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t)*j;
    
    nsc=100;%海明窗的长度,即每个窗口的长度
    nov=0;%重叠率
    nff= max(256,2^nextpow2(nsc));%N点采样长度
    %% matlab自带函数
    figure(1)
    spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
    title('spectrogram函数画图')
    [spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);

    变量:

    其中:

    spec_X矩阵行数为nff,列数是根据信号signal的长度和每个片段nsc分割决定的,共计col列。

    spec_f 为nff*1的列向量。

    spec_t为 1*ncol的行向量。

    此时自己实现的Matlab代码为:

      1 % 结合之前两个版本的stft,实现自己的周期图法,力求通俗易懂,代码分明。
      2 clear; clc; close all
      3 %% 信号输入基本参数
      4 Fs = 1000;            % Sampling frequency
      5 T = 1/Fs;             % Sampling period
      6 L = 1000;             % Length of signal
      7 t = (0:L-1)*T;        % Time vector
      8 S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t)*j;
      9 
     10 
     11 %传参
     12 signal = S;
     13 nsc    = 100;                     %每个窗口的长度,也即海明窗的长度
     14 nff    = max(256,2^nextpow2(nsc));%每个窗口进行fft的长度
     15 nov    = 0;                       %重叠率
     16 fs     = Fs;                      %采样率
     17 
     18 %% 使用fft实现周期图法
     19 %后续可封装为函数:
     20 % function [ S , W , T ] = mf_spectrogram...
     21 %     ( signal , nsc , nov , fs )
     22 % Input        
     23 %       signal      - Signal vector         输入信号(行向量)
     24 %       nsc         - Number SeCtion        每个小分段的信号长度
     25 %       nov         - Number OverLap        重叠点数
     26 %       fs          - Frequency of Sample   采样率
     27 % Output
     28 %       S           - A matrix that each colum is a FFT for time of nsc 
     29 %                    如果nfft为偶数,则S的行数为(nfft/2+1),如果nfft为奇数,则行数为(nfft+1)/2   
     30 %                    每一列是一个小窗口的FFT结果,因为matlab的FFT结果是对称的,只需要一半
     31 %       W           - A vector labeling frequency   频率轴
     32 %       T           - A vector labeling time        时间轴
     33 % Signal Preprocessing
     34 h = hamming(nsc, 'periodic');        % Hamming weight function  海明窗加权,窗口大小为nsc
     35 L = length(signal);                  % Length of Signal         信号长度
     36 nstep  = nsc-nov;                    % Number of STep per colum 每个窗口的步进
     37 ncol   = fix( (L-nsc)/nstep ) + 1;   % Number of CoLum          信号被分成了多少个片段
     38 nfft   = max(256,2^nextpow2(nsc));   % Number of Fast Fourier Transformation  每个窗口FFT长度
     39 nrow   = nfft/2+1;
     40 %
     41 %
     42 STFT1 = zeros(nfft,ncol);  %STFT1 is a matrix 初始化最终结果
     43 index  = 1;%当前片段第一个信号位置在原始信号中的索引
     44 for i=1:ncol
     45     %提取当前片段信号值,并用海明窗进行加权(均为长度为nsc的行向量)
     46     temp_S=signal(index:index+nsc-1).*h';
     47     %对长度为nsc的片段进行nfft点FFT变换
     48     temp_X=fft(temp_S,nfft);
     49     %将长度为nfft点(行向量)的fft变换转换为列向量,存储到矩阵的列向量
     50     STFT1(:,i)=temp_X';
     51     %将索引后移
     52     index=index + nstep;
     53 end
     54 
     55 
     56 %如果是为了和标准周期图法进行误差比较,则后续计算(abs,*2,log(+1e-6))不需要做,只需
     57 %plot(abs(spec_X)-abs(STFT1))比较即可
     58 
     59 STFT2 = abs(STFT1/nfft);
     60 %nfft是偶数,说明首尾是奈奎斯特点,只需对2:end-1的数据乘以2
     61 STFT2(2:end-1,:) = 2*STFT2(2:end-1,:);  % Add the value of the other half
     62 %STFT3 = 20*log10(STFT1);                % Turn sound pressure into dB level
     63 STFT3 = 20*log10(STFT2 + 1e-6);          % convert amplitude spectrum to dB (min = -120 dB)
     64 
     65 % Axis Generating
     66 fax =(0:(nfft-1)) * fs./nfft;                        % Frequency axis setting
     67 tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs;     % Time axis generating
     68  
     69  % Output
     70  W = fax;
     71  T = tax;
     72  S = STFT3;
     73  
     74 
     75 %% matlab自带函数
     76 figure();
     77 spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
     78 title('spectrogram函数画图')
     79 [spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
     80 %减法,看看差距
     81 figure()
     82 plot(abs(spec_X)-abs(STFT1))
     83 
     84 %% 画频谱图
     85 figure
     86 %利用了SIFT3
     87 surf(W, T,  S')
     88 shading interp
     89 axis tight
     90 box on
     91 view(0, 90)
     92 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
     93 xlabel('Frequency, Hz')
     94 ylabel('Time, s')
     95 title('Amplitude spectrogram of the signal')
     96 handl = colorbar;
     97 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
     98 ylabel(handl, 'Magnitude, dB')
     99 
    100 %跳频图
    101 figure();
    102 surf(T,W,S)
    103 shading interp
    104 axis tight
    105 box on
    106 view(0, 90)
    107 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
    108 ylabel('Frequency, Hz')
    109 xlabel('Time, s')
    110 title('Amplitude spectrogram of the signal')
    111 handl = colorbar;
    112 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
    113 ylabel(handl, 'Magnitude, dB')
    View Code

    重点关注8行、59行、66行、69行、82行

    其实此时只要牢牢记住结果集中行数为nfft就行了。

    2)当输入信号为实数时的代码:

    clear
    clc
    close all
    Fs = 1000; % Sampling frequency
    T = 1/Fs; % Sampling period
    L = 1000; % Length of signal
    t = (0:L-1)*T; % Time vector
    S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t);

    nsc=100;%海明窗的长度,即每个窗口的长度
    nov=0;%重叠率
    nff= max(256,2^nextpow2(nsc));%N点采样长度
    %% matlab自带函数
    figure(1)
    spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
    title('spectrogram函数画图')
    [spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);

     

    其中:

    spec_X矩阵行数为nff/2+1,列数是根据信号signal的长度和每个片段nsc分割决定的,共计col列。

    spec_f 为nff/2+1*1的列向量。

    spec_t为 1*ncol的行向量。

    主要是因为输入信号为实数时,fft变换的结果有一半是对称的,不必考虑。

    自己实现的Matlab代码为:

      1 % 结合之前两个版本的stft,实现自己的周期图法,力求通俗易懂,代码分明。
      2 clear; clc; close all
      3 %% 信号输入基本参数
      4 Fs = 1000;            % Sampling frequency
      5 T = 1/Fs;             % Sampling period
      6 L = 1000;             % Length of signal
      7 t = (0:L-1)*T;        % Time vector
      8 S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t);
      9 
     10 
     11 %传参
     12 signal = S;
     13 nsc    = 100;                     %每个窗口的长度,也即海明窗的长度
     14 nff    = max(256,2^nextpow2(nsc));%每个窗口进行fft的长度
     15 nov    = 0;                       %重叠率
     16 fs     = Fs;                      %采样率
     17 
     18 %% 使用fft实现周期图法
     19 %后续可封装为函数:
     20 % function [ S , W , T ] = mf_spectrogram...
     21 %     ( signal , nsc , nov , fs )
     22 % Input        
     23 %       signal      - Signal vector         输入信号(行向量)
     24 %       nsc         - Number SeCtion        每个小分段的信号长度
     25 %       nov         - Number OverLap        重叠点数
     26 %       fs          - Frequency of Sample   采样率
     27 % Output
     28 %       S           - A matrix that each colum is a FFT for time of nsc 
     29 %                    如果nfft为偶数,则S的行数为(nfft/2+1),如果nfft为奇数,则行数为(nfft+1)/2   
     30 %                    每一列是一个小窗口的FFT结果,因为matlab的FFT结果是对称的,只需要一半
     31 %       W           - A vector labeling frequency   频率轴
     32 %       T           - A vector labeling time        时间轴
     33 % Signal Preprocessing
     34 h = hamming(nsc, 'periodic');        % Hamming weight function  海明窗加权,窗口大小为nsc
     35 L = length(signal);                  % Length of Signal         信号长度
     36 nstep  = nsc-nov;                    % Number of STep per colum 每个窗口的步进
     37 ncol   = fix( (L-nsc)/nstep ) + 1;   % Number of CoLum          信号被分成了多少个片段
     38 nfft   = max(256,2^nextpow2(nsc));   % Number of Fast Fourier Transformation  每个窗口FFT长度
     39 nrow   = nfft/2+1;
     40 %
     41 %
     42 STFT1 = zeros(nfft,ncol);  %STFT1 is a matrix 初始化最终结果
     43 index  = 1;%当前片段第一个信号位置在原始信号中的索引
     44 for i=1:ncol
     45     %提取当前片段信号值,并用海明窗进行加权(均为长度为nsc的行向量)
     46     temp_S=signal(index:index+nsc-1).*h';
     47     %对长度为nsc的片段进行nfft点FFT变换
     48     temp_X=fft(temp_S,nfft);
     49     %将长度为nfft点(行向量)的fft变换转换为列向量,存储到矩阵的列向量
     50     STFT1(:,i)=temp_X';
     51     %将索引后移
     52     index=index + nstep;
     53 end
     54 
     55 % -----当输入信号为实数时,对其的操作(也就是说只考虑信号的前一半)
     56 % 由于实数FFT的对称性,提取STFT1的nrow行
     57 STFT_X = STFT1(1:nrow,:);             % Symmetric results of FFT
     58 
     59 %如果是为了和标准周期图法进行误差比较,则后续计算(abs,*2,log(+1e-6))不需要做,只需
     60 %plot(abs(spec_X)-abs(STFT_X))比较即可
     61 
     62 STFT2 = abs(STFT_X/nfft);
     63 %nfft是偶数,说明首尾是奈奎斯特点,只需对2:end-1的数据乘以2
     64 STFT2(2:end-1,:) = 2*STFT2(2:end-1,:);  % Add the value of the other half
     65 %STFT3 = 20*log10(STFT1);                % Turn sound pressure into dB level
     66 STFT3 = 20*log10(STFT2 + 1e-6);          % convert amplitude spectrum to dB (min = -120 dB)
     67 
     68 % Axis Generating
     69 fax =(0:(nrow-1)) * fs./nfft;                        % Frequency axis setting
     70 tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs;     % Time axis generating
     71  
     72  % Output
     73  W = fax;
     74  T = tax;
     75  S = STFT3;
     76  
     77 
     78 %% matlab自带函数
     79 figure();
     80 spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
     81 title('spectrogram函数画图')
     82 [spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
     83 %减法,看看差距
     84 figure()
     85 plot(abs(spec_X)-abs(STFT_X))
     86 
     87 %% 画频谱图
     88 figure
     89 %利用了SIFT3
     90 surf(W, T,  S')
     91 shading interp
     92 axis tight
     93 box on
     94 view(0, 90)
     95 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
     96 xlabel('Frequency, Hz')
     97 ylabel('Time, s')
     98 title('Amplitude spectrogram of the signal')
     99 handl = colorbar;
    100 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
    101 ylabel(handl, 'Magnitude, dB')
    102 
    103 %跳频图
    104 figure();
    105 surf(T,W,S)
    106 shading interp
    107 axis tight
    108 box on
    109 view(0, 90)
    110 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
    111 ylabel('Frequency, Hz')
    112 xlabel('Time, s')
    113 title('Amplitude spectrogram of the signal')
    114 handl = colorbar;
    115 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
    116 ylabel(handl, 'Magnitude, dB')
    View Code

    主要记到行数为nrow=nfft/2+1行。

    重点关注代码8行,57行,62行,85行。

    二、分析思路

    1)自己在使用fft实现时,通过一个for循环,对每个nsc长度的片段进行加窗、nfft点的fft变换。

    2)

    1. 当输入信号为实数时,由于fft的对称性,从结果SIFT1仅取fft结果行数的一半放到结果矩阵SIFT_X中去。

    2. 当输入信号为复数时,该矩阵的行数仍然为nff行,因此结果集即为SIFT1。

    3. 此时每一列即是一个小片段的fft变换的结果,该列结果是复数。 可以和Matlab API得到的结果进行差值比较,发现结果完全一样。

    这3处分析的不同,也正是两个自己实现代码的修改处,可参考上述代码。

    3)后续由于考虑到绘图、幅值、频谱的影响,因此又在复数矩阵的结果上进行了一些运算,因此在后续封装为函数时根据需要进行改进。

     

    转载于:https://www.cnblogs.com/shuqingstudy/p/10155980.html

    展开全文
  • Matlab中,做短时傅里叶变换需要使用函数spectrogram,而在Matlab2019中,引入了一个新的函数stft,下面我们就来看下这两个函数都如何使用。 短时傅里叶变换的基本原理就是将数据分段加窗,做fft,在分段时会有...
  •   在Matlab中,做短时傅里叶变换需要使用函数spectrogram,而在Matlab2019中,引入了一个新的函数stft,下面我们就来看下这两个函数都如何使用。  短时傅里叶变换的基本原理就是将数据分段加窗,做fft,在分段时...
  • Matlab短时傅里叶变换 spectrogram和stft的用法

    千次阅读 多人点赞 2020-05-25 19:22:29
      在Matlab中,做短时傅里叶变换需要使用函数spectrogram,而在Matlab2019中,引入了一个新的函数stft,下面我们就来看下这两个函数都如何使用。   短时傅里叶变换的基本原理就是将数据分段加窗,做fft,在分段...
  • matlab时频分析之短时傅里叶变换 spectrogram

    万次阅读 多人点赞 2019-03-22 15:50:25
    matlab时频分析之短时傅里叶变换 spectrogram 短时傅里叶变换常用于缓慢时变信号的频谱分析,可以观察沿时间变化的频谱信号。 其优点如下图所示,弥补了频谱分析中不能观察时间的缺点,也弥补了时域分析不能获取频率...
  • Matlab短时逆傅里叶变换函数 本节内容部分来自于Matlab参考文档 Matlab中短时逆傅里叶函数是signal processing toolbox工具箱中的istft函数,调用格式x = istft(s,fs,Name,Value); 若想使istft后的函数完美复现原...
  •   在Matlab中,做短时傅里叶变换需要使用函数spectrogram,而在Matlab2019中,引入了一个新的函数stft,下面我们就来看下这两个函数都如何使用。  短时傅里叶变换的基本原理就是将数据分段加窗,做fft,在分段时...
  • 摘要:一直以来都是用MATLAB做信号处理,得到预处理的特征后再用Python...这里将MATLAB中的FFT、STFT、加窗以及带通滤波通过Python接口实现,防止以后MATLAB用不了了,一定程度上也提高了效率,不用两个软件换来换去。
  • % 线性调频信号% clear all k=50; fs=250; n=0:1000-1; x=exp(j*pi*k(n/fs^2; subplot(221; plot(n/fs,real(x; title'线性调频信号;... % 线性调频信号频谱 % P=fft(x; subplot(222; plot(n/length(n*fs,abs(P; title
  • Matlab绘制短时傅里叶变换的频谱图和时间-频率-幅值三维图 function [t,frequency,f_spectrum]=fft_s(y,windowlength,Fs) % 输入 : % y-输入信号 % windowlength-窗长度 0-1的系数,比如windowlength-0.5;采样...
  • 数字信号处理(3)- 短时频域分析

    千次阅读 2021-01-23 19:56:45
    当N是2的整数倍时,这个离散短时傅里叶变换可以使用FFT来计算。 MATLAB程序 MATLAB程序演示信号分帧、加窗、求离散短时傅里叶变换,并最终使用三维图展示结果。 其中打开的test.wav文件是一个8kHz采样率的音频文
  • MATLAB不同时频信号处理方法介绍及效果对比

    千次阅读 多人点赞 2018-10-25 17:29:40
    目录MATLAB不同时频信号处理方法介绍及效果对比信号的频域分析Fourier(FFT)变换基本原理MATLAB 中信号FFT的基本处理方法MATLAB 中不同参数设置对FFT结果的影响信号的功率谱分析的基本原理MATLAB中信号功率谱分析的...
  • 该法可对语音信号(或其他类似的似平稳信号)进行基于FFT短时频谱分析,频谱图的伪彩色映射及显示。频谱图的类型有宽带和窄带两种,伪彩色显示的映射可按需求设置。文中介绍了基于FFT频谱分析与频谱图显示的原理,...
  • STFT短时傅里叶变换的实现

    万次阅读 2017-04-24 18:53:59
    本科毕设做雷达微动的时频分析,最常用的就是STFT了。matlab已经有了时频分析工具箱,调用里面的tfrstft函数即可轻松实现。今天看了源代码感觉写的比较复杂,于是结合微动,尝试自己实现一下。  STFT是在...
  • 1、FFT,分析基波的参与时间对傅里叶变换的影响 基波为4个余弦波,用matlab程序控制其参与输入信号的...注意:FFT后在Vector Scope显示每帧的频率图,应设置采样频率,否则,频率计算可能有误! 比如,Buffer设...
  • 分帧是为了将无限长的语音信号,分成一段一段的,因为语音信号具有短时平稳性,方便处理,加窗是为了使分帧后的语音信号更加平稳。窗函数主要有矩形窗和汉明窗。加窗主要是为了使时域信号似乎更好地满足FFT处理的...
  • 写出一段matlab程序,实现下述功能: 1,读取作业一的录音文件 2,...对该帧进行短时幅度和FFT频谱变换,并画图显示。 上述三个图像,用subplot命令, 画在一张图里</p>
  • 叹息......我没有给出Matlab ...首先,您要使用FFT将时域信号转换为短时频域(STFT):这可以通过matlab specgram 来完成:[data,fs,bit] = wavread('test.wav'); %//Read in WAV file%//Note if you have dual ch...
  • :) FFT频谱分析 短时傅里叶变换 经典功率谱估计(weltch法和Thomson多窗估计法) 现代功率谱估计AR模型法 1.1.jpg 1.1tu.jpg 1.2.jpg...
  • 语音识别的MATLAB实现

    热门讨论 2009-03-03 21:39:18
    在计算短时能量之前应用该滤波器,还可以起到消除直流漂移、抑制随机噪声和提升清音部分能量的效果。 C 分帧处理 在计算各个系数之前要先将语音信号作分帧处理。语音信号是瞬时变化的,但在10~20ms内是相对稳定...
  • 该函数的功能是将时域的语音信号变为短时频域信号; 主要包括:预加重,分帧,加窗,FFT求能量值 designAuditoryFilterBank designAuditoryFilterBank 该函数的功能是将频域的语音信号通过滤波器组进行滤波,如常见...
  • 语音信号为一维短时平稳信号,带宽为3~4kHz,分析时可按10~20ms相对稳定的时间段分段进行,对每一段信号可直接采用一维FFT算法分析其频谱,或描述其对应的语谱图,即横轴为时间、纵轴为频率和图像的灰度表示相应...
  • 本文不涉及MFCC的理论,所以读此文前请对MFCC以及相关语音信号处理有初步认识。本文重点在于代码实现的分析。 先对MFCC有个初步认识。 MFCCs(Mel Frequency... 2)对每一个短时分析窗,通过FFT得到对应的频谱; 3)...
  • 2.分帧计算语音信号倒谱,并倒谱作FFT并加短时窗。 3.取大于25以上的样值,进行IFFT,得到基音周期的倒谱。 3.运用Levinson-Durbin计算一帧语音信号线形预测系数 4.进行逆滤波处理 5.对逆滤波后的信号进行倒谱...

空空如也

空空如也

1 2
收藏数 36
精华内容 14
热门标签
关键字:

matlab短时fft

matlab 订阅