# 数字信号处理公式变程序（五）——仿matlab的spectrogram函数（STFT）

STFT开始！

============================================================

1.概念和特点

STFT（short-time Fourier transform，短时傅里叶变换）是和傅里叶变换相关的一种数学变换，用以确定时变信号其局部区域正弦波的频率与相位。

2.原理简介

3.设计思路

(1)窗函数选择hamming窗，最大DFT点数不大于256；

(2)用户输入（传值）：signal, window, overlap, N, fs等；

(3)根据窗的大小，将signal拆分，并与窗函数相乘；

(4)对每个signal片段进行N点FFT，并求出能量谱密度；

(5)调用绘图方法，把能量谱密度（功率谱密度）用不同的颜色表示出来绘图。

(1)各种窗表达式

①海明窗(hamming)：w(n)=0.54-0.46*cos(n/N);

②汉宁窗(hanning)：w(n)=0.5*(1-cos(n/N));

③矩形窗(Rectangular)：w(n)=1.0;

④三角窗(Triangle)：w(n)=TRI(2n/N);

⑤布莱克曼窗(Blackman,三阶升余弦窗)：w(n)=0.42-0.5*cos(n/N)+0.08*cos(2n/N);

⑥布莱克曼-哈里斯窗(BlackmanHarris)：w(n)=0.35875-0.48829*cos(n/N)+0.14128*cos(2n/N)-0.01168*cos(3n/N);

(2)信号与窗的相乘

(3)绘图用到的颜色

①在matlab中colorbar可以调出颜色条图，如colorbar,colorbar('North')等，North表示在top出现颜色条，另外还有East, West, South等（详见matlab帮助help colorbar）；

②在matlab中colormap可以调出颜色条图对应的RGB值数组；

③复制一下RGB对应的数组如下：

#pragma mark - 颜色数组

//colorbar,colormap，定义冷暖色条的RGB值
static float g_colorBarR[64] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.06250,0.1250,0.18750,0.2500,0.31250,0.3750,0.43750,0.500,0.56250,0.6250,0.68750,0.7500,0.81250,0.8750,0.93750,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.93750,0.8750,0.81250,0.7500,0.68750,0.6250,0.56250,0.500};
static float g_colorBarG[64] = {0,0,0,0,0,0,0,0,0.0625,0.125,0.1875,0.250,0.3125,0.375,0.4375,0.500,0.5625,0.625,0.6875,0.750,0.8125,0.875,0.9375,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.9375,0.875,0.8125,0.750,0.6875,0.625,0.5625,0.500,0.4375,0.375,0.3125,0.250,0.1875,0.125,0.0625,0,0,0,0,0,0,0,0,0};
static float g_colorBarB[64] = {0.5625,0.625,0.6875,0.750,0.8125,0.875,0.9375,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.9375,0.875,0.8125,0.750,0.6875,0.625,0.5625,0.500,0.4375,0.375,0.3125,0.250,0.1875,0.125,0.0625,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
(4)能量谱密度（功率谱密度PSD）

①设一个能量信号 s(t)

②在MATLAB中计算过程是

1.spectrogram方法

/**
*  @brief   模仿MATLAB的spectrogram，实现STFT功能
*  @param   signalVector 原始信号序列
*  @param   N            信号的长度
*  @param   win          窗长度
*  @param   noverlap     窗重叠的点数
*  @param   nfft         FFT/DFT的点数
*  @param   fs           抽样频率
*  @param   drawRect     绘图的区域rect
*  @param   context      绘图的上下文
*  @return  是否计算成功的标志
*/
+ (Boolean)spectrogram:(float *)signalVector andLength:(int)N andWindowLength:(int)win andNoverlap:(int)noverlap andDFTCount:(int)nfft andFs:(int)fs andDrawRect:(CGRect)drawRect andContext:(CGContextRef)context
{
Boolean isSpecOK = false;
//hamming窗计算
float hammingWindow[win], windowPowValue = 0.0;
for(int i = 0; i < win; i++)
{
hammingWindow[i] = 0.54 - 0.46 * cosf(2*M_PI*i/(win-1));
windowPowValue += powf(hammingWindow[i], 2);
}

//计算短时傅里叶变换信号数组的行与列数，行数=时间点数，列数=窗长度
int row, column, halfNfft;
row = (N - noverlap)/(win - noverlap);
column = win;
halfNfft = nfft/2+1;

//计算t数组
float timeVector[row];
for(int i = 0; i < row; i++)
{
timeVector[i] = ((float)i) /((float)(fs*(win/2+1+(win-noverlap)*i)));
}

//将signal拆分
float signalXY[row][column];
for(int i = 0; i < row; i++)
{
for(int j = 0; j < column; j++)
{
signalXY[i][j] = signalVector[i*(win - noverlap) + j];
signalXY[i][j] *= hammingWindow[j];
}
}

//对拆分后的信号进行FFT
float fftReal[nfft], fftImg[nfft], Sxx[row][halfNfft], Pxx[row][halfNfft], logPxx[row][halfNfft], freq[nfft], freqStep, pxxMax, pxxMin;

freqStep = ((float)fs)/((float)nfft);
freq[0] = 0.0;
for(int i = 1; i < nfft; i++)
{
freq[i] = freqStep + freq[i-1];
}

//求Sxx
for(int i = 0; i < row; i++)
{
[FFT fft:signalXY[i] andOriginalLength:column andFFTCount:nfft andFFTReal:fftReal andFFTYImg:fftImg];
for(int j = 0; j < halfNfft; j++)
{
Sxx[i][j] = (fftReal[j]*fftReal[j] + fftImg[j]*fftImg[j])/windowPowValue;
}
}

//求Pxx
float fsFloat = (float)fs;
Pxx[0][0] = Sxx[0][0]/fsFloat;
logPxx[0][0] = 10*log10f(fabsf(Pxx[0][0]));
pxxMax = logPxx[0][0];
pxxMin = logPxx[0][0];

for(int i = 1; i < row; i++)
{
Pxx[i][0] = Sxx[i][0]/fsFloat;
logPxx[i][0] = 10*log10f(fabsf(Pxx[i][0]));

if(logPxx[i][0] > pxxMax)
pxxMax = logPxx[i][0];
if(logPxx[i][0] < pxxMin)
pxxMin = logPxx[i][0];
}
if(nfft%2)//奇数
{
for(int i = 0; i < row; i++)
{
for(int j = 1; j < halfNfft; j++)
{
Pxx[i][j] = Sxx[i][j]*2.0/((float)fs);
logPxx[i][j] = 10*log10f(fabsf(Pxx[i][j]));

if(logPxx[i][j] > pxxMax)
pxxMax = logPxx[i][j];
if(logPxx[i][j] < pxxMin)
pxxMin = logPxx[i][j];
}
}
}
else//偶数
{
for(int i = 0; i < row; i++)
{
Pxx[i][halfNfft-1] = Sxx[i][halfNfft-1]/((float)fs);
logPxx[i][halfNfft-1] = 10*log10f(fabsf(Pxx[i][halfNfft-1]));

if(logPxx[i][halfNfft-1] > pxxMax)
pxxMax = logPxx[i][halfNfft-1];
if(logPxx[i][halfNfft-1] < pxxMin)
pxxMin = logPxx[i][halfNfft-1];
}
for(int i = 0; i < row; i++)
{
for(int j = 1; j < halfNfft-1; j++)
{
Pxx[i][j] = Sxx[i][j]*2.0/((float)fs);
logPxx[i][j] = 10*log10f(fabsf(Pxx[i][j]));

if(logPxx[i][j] > pxxMax)
pxxMax = logPxx[i][j];
if(logPxx[i][j] < pxxMin)
pxxMin = logPxx[i][j];
}
}
}

//绘图
NSMutableArray *dataArray = [[NSMutableArray alloc] init];
for(int i = 0; i < row; i++)
{
for(int j = 0; j < halfNfft; j++)
{
NSNumber *numberMax = [NSNumber numberWithFloat:logPxx[i][j]];
}
}
NSArray *colorArray = [[NSArray alloc] initWithArray:[Drawing getColorArrayByValue:dataArray andLengthX:row andLengthY:halfNfft andMax:pxxMax andMin:pxxMin]];

CGPoint p1, p2;
p1 = drawRect.origin;
p2.x = p1.x + drawRect.size.width;
p2.y = p1.y + drawRect.size.height;
Drawing *drawing = [[Drawing alloc] init];
[drawing setEnableDrawFieldLeftTopPoint:p1 andRightBottomPoint:p2 andContext:context];

int nfftAll = [FFT nextNumOfPow2:N];
float fftRealTest[nfftAll], fftImgTest[nfftAll];

[FFT fft:signalVector andOriginalLength:N andFFTCount:nfftAll andFFTReal:fftRealTest andFFTYImg:fftImgTest];

float t[N];
for(int i = 0; i < N; i++)
{
printf("%f + %fi   ", fftRealTest[i], fftImgTest[i]);
t[i] = i;
fftRealTest[i] = sqrtf((fftImgTest[i]*fftImgTest[i] + fftRealTest[i]*fftRealTest[i]))/nfftAll;
}

[drawing selectSubDrawFieldWithRow:3 andColumn:1 andNumber:1];
[drawing drawLineWithXVector:t andYVector:signalVector andLength:N andLineColor:[UIColor blueColor]];

[drawing selectSubDrawFieldWithRow:3 andColumn:1 andNumber:2];
[drawing drawLineWithXVector:t andYVector:fftRealTest andLength:nfftAll/2 andLineColor:[UIColor blueColor]];

[drawing selectSubDrawFieldWithRow:3 andColumn:1 andNumber:3];
[drawing drawColorPointsWithLengthX:row andLengthY:halfNfft andColorArray:colorArray];

return isSpecOK;
}
2.绘图方法

//一次性画多个颜色点
- (void)drawColorPointsWithLengthX:(int)lx andLengthY:(int)ly andColorArray:(NSArray *)colors
{
if([colors count] != lx*ly)
{
NSLog(@"错误，长度不匹配！");
exit(1);
}

float point_x, point_y;
UIColor *pointColor = [colors objectAtIndex:0];
float divX, divY, widthX, widthY;

divX = (cornerPoints[1].x - cornerPoints[0].x - 0.6)/((float)lx);
divY = (cornerPoints[2].y - cornerPoints[0].y - 1.0)/((float)ly);
widthX = ceilf(divX);
widthY = ceilf(divY);

for(int i = 0; i < lx; i++)
{
point_x = cornerPoints[0].x + 0.6 + divX*i;
int countx = i*ly;
for(int j = 0; j < ly; j++)
{
point_y = cornerPoints[2].y - 2.5 - divY*j;
pointColor = [colors objectAtIndex:countx+j];

CGContextSetFillColorWithColor(drawContext, [pointColor CGColor]);
CGContextFillEllipseInRect(drawContext, CGRectMake(point_x, point_y, widthX, widthY));
//CGContextFillPath(drawContext);  //填充颜色
}
}
}

1.先看看STFT观察跳频

signal1(freq=100Hz,Amp=2V)

signal2(freq=150Hz,Amp=2V)

signal3(freq=300Hz,Amp=1.5V)

①利用MATLAB测试

%% 变频信号的fft观察
close all
clear,clc
%% 定义一个变频信号的参数
PIx2 = 2 * pi;
Fs = 1500;
T = 1/Fs;
LengthOfSignal = 3000;
NFFT = 2^nextpow2(LengthOfSignal); %要计算的点数
% NFFT = 128;
t = (0:LengthOfSignal-1)*T;
amp1 = 2;
amp2 = 2;
amp3 = 1.5;
offset = 2;
freq1 = 100;
freq2 = 150;
freq3 = 300;
signal = zeros(1,length(t));

%% 定义信号
for temp = 1:LengthOfSignal
if(temp <= LengthOfSignal/4)
signal(temp) =offset + amp1 * sin(PIx2 * freq1 * t(temp));
elseif(temp <= LengthOfSignal/2)
signal(temp) =offset -1*amp2 * sin(PIx2 * freq2 * t(temp));
elseif(temp <= 3*LengthOfSignal/4)
signal(temp) =offset -1*amp3 * sin(PIx2 * freq3 * t(temp));
else
signal(temp) =offset + amp1 * sin(PIx2 * freq1 * t(temp));
end
end

%% 绘制图形
subplot(311);
plot(t, signal);
grid on
title('signal of different frequency');
xlabel('time');
ylabel('amp');

%% fft运算
fMax = NFFT/2 + 1;
signalFFT = abs(fft(signal,NFFT));
% signalFFTShow = 2 * abs(fft(signal(1:fMax),NFFT)/LengthOfSignal);
signalFFTShow = 2 * signalFFT / LengthOfSignal;
signalFFTShow(1) = signalFFTShow(1)/2;
f = Fs/2*linspace(0,1,fMax);
subplot(312);
plot(f,signalFFTShow(1:fMax));
grid on
title('fft signal');
xlabel('frequency');
ylabel('amp');

%% surf测试三维图(不合理x(j),y(i),z(i,j)对应)

subplot(313)
spectrogram(signal,128,120,128,1.5e3,'yaxis'); %时频域图


③结果简单总结,从上面两幅图可以看出，

A.matlab运行的结果含有直流分量绘图（iOS去掉了）。

B.另外颜色有偏差，原因可能是功率谱密度与冷暖色条对应的方式不同（我没深究。。。）。

2.再看看STFT观察对语音的分析

①看看matlab程序及结果

%% spectrogram分析chirp声音
close all
clear,clc

%% 基本运算部分
signalOfChirp = chirpData';
signalOfChirp_5 = signalOfChirp/5;
lengthOfSignal = length(signalOfChirp);
Fs = lengthOfSignal;  %采样频率
T = 1/Fs;
t = (0:lengthOfSignal-1)*T;

% 绘原始波形图
figure;
subplot(311);
plot(t, signalOfChirp);
grid on
title('Original signal(chirpData)');
xlabel('Time');
ylabel('Amp');

%% 计算信号的fft
NFFT = 2^nextpow2(lengthOfSignal); %求最接近2的幂的数作为FFT计算的点数
fMax = NFFT/2 + 1;                 %频域显示最大的频率
f = Fs/2*linspace(0,1,fMax);       %确定f数组

signalFFT = abs(fft(signalOfChirp, NFFT)/lengthOfSignal); %求出fft并取模
signalFFT_5 = 5 * abs(fft(signalOfChirp_5, NFFT)/lengthOfSignal);
%绘制fft波形图
subplot(312);
plot(f,signalFFT(1:fMax));
grid on
title('FFT wave of signal');
xlabel('Frequency');
ylabel('Amp');

%% 绘制f-t图（STFT）
subplot(313);

spectrogram(signalOfChirp,128,120,128,Fs,'yaxis'); %时频域图

title('f-t');

②看看ios的测试结果

③结论

3.比一比Matlab和iOS的效率如何

    NSDate *start1 = [NSDate date];

// ...需要测试的代码

NSDate *end1 = [NSDate date];
NSLog(@"FFT time %f ", [end1 timeIntervalSinceDate:start1]);

============================================

OK，STFT也结束了~~

• 本文已收录于以下专栏：

举报原因： 您举报文章：数字信号处理公式变程序（五）——仿matlab的spectrogram函数（STFT） 色情 政治 抄袭 广告 招聘 骂人 其他 (最多只允许输入30个字)