Matlab周期图法使用FFT实现
2024-05-10 15:25:07
参考文章：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 dBSTFT1 = abs(STFT_X/nfft);STFT1(2:end-1,:) = 2*STFT1(2:end-1,:); % Add the value of the other halfSTFT3 = 20*log10(STFT1); % Turn sound pressure into dB level% Axis Generatingfax =(0:(nfft/2)) * fs/nfft; % Frequency axis settingtax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs; % Time axis generating[ffax,ttax] = meshgrid(tax,fax); % Generating grid of figure% OutputW = 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 dBs1 = abs(s/nff);s2 = s1(1:nff/2+1,:); % Symmetric results of FFTs2(2:end-1,:) = 2*s2(2:end-1,:); % Add the value of the other halfs3 = 20*log10(s2); % Turn sound pressure into dB levelmy_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周期图法使用FFT实现相关推荐
- matlab周期图法,周期图谱估计及其改进方法的Matlab仿真
1引言对信号与系统的研究处理主要有时域.频域两种方法.对确定性信号,可以进行傅里叶变换,从而进行频域分析,但对随机信号,由于其傅里叶变换不存在,通常是求其功率谱来进行频谱分析,因为功率谱反映了随机信号 ...
- matlab中拟合函数中的gian值,如何在Matlab中优化基本周期图法对随机信号进行的功率谱估计...
首都师范大学学报(自然科学版)第27卷 第5期2006年10月 Journal of Capital N ormal University (Natural Science Edition ) V o ...
- matlab fft实现dft,matlab实现dft和fft
对任意长度的序列进行傅里叶变换 DFT 与 FFT 的运算时间比较 设计要求 利用 Matlab 或者 C 语言设计 DFT 和 FFT 程序,比较两种频谱分析方法的 计算速度,并与...... DF ...
- matlab晶体能带,matlab平面波展开法的二维光子晶体能带研究+程序
摘 要 :二维光子晶体可以作为对光子传输控制的新型材料.本文主要通过平面波展开法对二维光子晶体进行数值计算及其性质分析.首先我们介绍了二维光子晶体的基础概念.结构.介电性能等特性.然后基于麦克斯韦方 ...
- 【聚类算法】基于matlab划分法k-means聚类算法【含Matlab源码 1941期】
一.获取代码方式 获取代码方式1: 完整代码已上传我的资源:[聚类算法]基于matlab划分法k-means聚类算法[含Matlab源码 1941期] 点击上面蓝色字体,直接付费下载,即可. 获取代码 ...
- matlab负序,Matlab对称分量法计算正序负序零序
Matlab对称分量法计算正序负序零序 (2012-11-30 09:15:38) 标签: 分量 电压 杂谈 clc clear all close all prompt{1}='A相电压:'; pr ...
- MATLAB 欧几里得算法
MATLAB 欧几里得算法以及扩展欧几里得算法 前言 一.Euclid算法 二.Code 1.引入库 2.读入数据 前言 提示:在线性代数或高等代数中以熟悉求两个多项式最大公因子 一.Euclid算法 ...
- newmark法 matlab,newmark法程序newmark法程序.doc
newmark法程序newmark法程序 用matlab编程法 一.法原理 Newmark-?法是一种逐步积分的方法,避免了任何叠加的应用,能很好的适应非线性的反应分析. Newmark-?法假定: ...
- matlab 限幅,你好，不好意思打扰了，我在论坛里看到你回了一个用MATLAB限幅法编程降低PAPR，可以请教下您吗...
你好,不好意思打扰了,我在论坛里看到你回了一个用MATLAB限幅法编程降低PAPR,可以请教下您吗 更新时间:2019-04-15 19:12:26 107次访问 % papr of the ofdm ...
最新文章
- 居然以反战的名义，利用开源工具包“投毒”！
- 利用lrz、lsz工具在linux与windows之间传输文件
- charles 如何抓取web界面的包
- 面向对象编程（OOP）特性 类和对象
- 利用StringUtils可以避免空指针问题
- abp执行mysql语句_在ABP模板工程中使用MySql
- 【白皮书分享】2022新职业教育洞察白皮书：“职”成机遇，“育”见未来.pdf...
- Python实现超像素分割
- 计算机基础到底是哪些基础？为什么很重要！
- 编码器SRT协议三种模式(listener, caller, rendezvous)简介
- 你若运行，便是晴天！
- php中=%3e -%3e的区别,华为nova3e与nova2s买哪个好?华为nova2s和nova3e区别对比详细评测...
- 邓应海：焦点转向美国数据，央行会议！最新黄金走势分析
- 为什么Python是数据可视化编程的最佳选择？
- 3.神经网络-深度学习入门
- THREE.OrbitControls参数控制解释
- C#入门4——计算自由落体运动
- 基于JavaSwing开发画图软件 课程设计 小作业 大作业源码
- 可nbsp;爱nbsp;女nbsp;人
- python使用二分查找法查找某个数是否在列表内