常见的非矩形窗函数似乎都是对称的。有没有人想在 FFT 之前使用非对称窗函数的情况？（假设 FFT 孔径一侧的数据被认为比另一侧的数据更重要，或者噪音更小，等等）
如果是这样，研究了哪些类型的非对称窗函数，与（更有损？）偏移对称窗相比，它们将如何影响频率响应？
具有非对称窗口的 FFT？
信息处理
fft
傅里叶变换
窗函数
2021-12-29 01:39:02
3个回答
我将使用“窗口函数”的简写窗口。
对于音频，任何产生类似于预振铃或预回声的处理都会听起来像低比特率的 mp3 一样松散。当瞬态或脉冲的局部能量在时间上向后传播时会发生这种情况，例如通过修改重叠变换中的频谱数据，例如重叠修正离散余弦变换 (MDCT)。在这样的处理中，音频通过重叠的分析窗口进行窗口化，在频域中进行转换和处理（如将数据压缩到更小的比特率），再次使用合成窗口进行窗口化，然后将它们加在一起。分析和综合窗口的乘积必须使得重叠窗口总和为一。
传统上使用的窗函数是对称的，它们的宽度是频率选择性（长窗）和时域伪影避免（短窗）之间的折衷。窗口越宽，处理可以越早地传播信号。最近的解决方案是使用非对称窗口。使用的两个窗口可以是彼此的镜像。分析窗口从峰值快速下降到零，因此不会提前“检测到”脉冲，而合成窗口从零快速上升到峰值，因此任何处理的效果都不会在时间上向后传播太多。这样做的另一个优点是低延迟。非对称窗口可以具有良好的频率选择性，可以在音频压缩中替代可变大小的对称窗口，是一种万灵药。看M. Schnell、M. Schmidt、M. Jander、T. Albert、R. Geiger、V. Ruoppila、P. Ekstrand、M. Lutzky、B. Grill，“MPEG-4 增强型低延迟 AAC - 高音的新标准质量通信”，第 125 届 AES 大会，美国加利福尼亚州旧金山，预印本 7503，2008 年 10 月，以及另一篇会议论文，其中也展示了其窗口的傅立叶变换幅度： Schnell, M. 等。2007.增强的 MPEG-4 低延迟 AAC - 低比特率高质量通信。在第 122 届 AES 公约中。
图 1. 在重叠分析-处理-综合中使用不对称窗口的图示。分析窗口（蓝色）和合成窗口（黄橙色）的乘积（黑色虚线）与前一帧的窗口（灰色虚线）相加为一。使用 MDCT 时需要进一步的约束来保证完美的重建。
可以使用离散傅里叶变换（DFT、FFT）代替 MDCT，但在这种情况下会给出冗余的光谱数据。与 DFT 相比，MDCT 仅提供一半的光谱数据，同时如果选择合适的窗口仍能实现完美的重建。
这是我自己的不对称窗口设计（图 2），适用于使用 DFT 进行重叠分析-处理-合成，但不适用于无法提供完美重建的 MDCT。该窗口试图最小化均方时间和频率带宽的乘积（类似于受限的高斯窗口），同时保留一些潜在有用的时域属性：非负、单峰，峰值在“时间零”处，分析和合成在此附近窗口是彼此的镜像，函数和一阶导数连续，当窗口函数的平方被解释为非归一化概率密度函数时为零均值。该窗口使用差分进化进行了优化。
图 2. 左：适用于重叠分析-处理-再合成的不对称分析窗口及其时间反转对应合成窗口。右：余弦窗口，具有与非对称窗口相同的延迟
图 3. 图 2 的余弦窗口（蓝色）和非对称窗口（橙色）的傅里叶变换幅度。非对称窗口显示出更好的频率选择性。
这是绘图和不对称窗口的 Octave 源代码。绘图代码来自Wikimedia Commons。在 Linux 上，我建议先安装gnuplot, epstool, pstoedit,transfig并librsvg2-bin使用display.
pkg load signal
graphics_toolkit gnuplot
set (0, "defaultaxesfontname", "sans-serif")
set (0, "defaultaxesfontsize", 12)
set (0, "defaultaxeslinewidth", 1)
function plotWindow (w, wname, wfilename = "", wspecifier = "", wfilespecifier = "")
M = 32; % Fourier transform size as multiple of window length
Q = 512; % Number of samples in time domain plot
P = 40; % Maximum bin index drawn
dr = 130; % Maximum attenuation (dB) drawn in frequency domain plot
N = length(w);
B = N*sum(w.^2)/sum(w)^2 % noise bandwidth (bins)
k = [0 : 1/Q : 1];
w2 = interp1 ([0 : 1/(N-1) : 1], w, k);
if (M/N < Q)
Q = M/N;
endif
figure('position', [1 1 1200 600])
subplot(1,2,1)
area(k,w2,'FaceColor', [0 0.4 0.6], 'edgecolor', [0 0 0], 'linewidth', 1)
if (min(w) >= -0.01)
ylim([0 1.05])
set(gca,'YTick', [0 : 0.1 : 1])
else
ylim([-1 5])
set(gca,'YTick', [-1 : 1 : 5])
endif
ylabel('amplitude')
set(gca,'XTick', [0 : 1/8 : 1])
set(gca,'XTickLabel',[' 0'; ' '; ' '; ' '; ' '; ' '; ' '; ' '; 'N-1'])
grid('on')
set(gca,'gridlinestyle','-')
xlabel('samples')
if (strcmp (wspecifier, ""))
title(cstrcat(wname,' window'), 'interpreter', 'none')
else
title(cstrcat(wname,' window (', wspecifier, ')'), 'interpreter', 'none')
endif
set(gca,'Position',[0.094 0.17 0.38 0.71])
H = abs(fft([w zeros(1,(M-1)*N)]));
H = fftshift(H);
H = H/max(H);
H = 20*log10(H);
H = max(-dr,H);
k = ([1:M*N]-1-M*N/2)/M;
k2 = [-P : 1/M : P];
H2 = interp1 (k, H, k2);
subplot(1,2,2)
set(gca,'FontSize',28)
h = stem(k2,H2,'-');
set(h,'BaseValue',-dr)
xlim([-P P])
ylim([-dr 6])
set(gca,'YTick', [0 : -10 : -dr])
set(findobj('Type','line'),'Marker','none','Color',[0.8710 0.49 0])
grid('on')
set(findobj('Type','gridline'),'Color',[.871 .49 0])
set(gca,'gridlinestyle','-')
ylabel('decibels')
xlabel('bins')
title('Fourier transform')
set(gca,'Position',[0.595 0.17 0.385 0.71])
if (strcmp (wfilename, ""))
wfilename = wname;
endif
if (strcmp (wfilespecifier, ""))
wfilespecifier = wspecifier;
endif
if (strcmp (wfilespecifier, ""))
savetoname = cstrcat('Window function and frequency response - ', wfilename, '.svg');
else
savetoname = cstrcat('Window function and frequency response - ', wfilename, ' (', wfilespecifier, ').svg');
endif
print(savetoname, '-dsvg', '-S1200,600')
close
endfunction
N=2^17; % Window length, B is equal for Triangular and Bartlett from 2^17
k=0:N-1;
w = -cos(2*pi*k/(N-1));
w .*= w > 0;
plotWindow(w, "Cosine")
freqData = [0.66697133904805994131, -0.20556692772918355727, 0.49267389481655493588, -0.25062332863369246594, -0.42388422228212319087, 0.42317609537724842905, -0.03930334287740060856, -0.11936153294075849129, 0.30201210285940127687, -0.15541616804857899536, -0.16208119255594669039, 0.12843871362286504723, -0.04470810646117385351, -0.00521885027256757845, 0.07185811583185619522, -0.02835116723496184862, -0.01393644785822748498, 0.00780746224568363342, -0.00748496824751256583, 0.00119325723511989282, 0.00194602547595042175];
freqData(1) /= 2;
scale = freqData(1) + sum(freqData.*not(mod(1:length(freqData), 2)));
freqData /= scale;
w = freqData(1)*ones(1, N);
for bin = 1:(length(freqData)/2)
w += freqData(bin*2)*cos(2*pi*bin*((1:N)-1)/N);
w += freqData(bin*2+1)*sin(2*pi*bin*((1:N)-1)/N);
endfor
w(N/4+1:N/2+1) = 0;
w(N/8+2:N/4) = (1 - w(N/8:-1:2).*w(7*N/8+2:N))./w(7*N/8:-1:6*N/8+2);
w = shift(w, -N/2);
plotWindow(w, "Asymmetrical");
您可能只想使用窗口的每一秒样本，因为它从零开始和结束。以下 C++ 代码会为您执行此操作，因此您不会得到任何零样本，除非在窗口的四分之一处处处为零。对于分析窗口，这是第一季度，对于综合窗口，这是最后一个季度。分析窗口的后半部分应与合成窗口的前半部分对齐，以计算它们的乘积。该代码还测试了窗口的平均值（作为概率密度函数），并展示了重叠重建的平坦度。
#include <stdio.h>
#include <math.h>
int main() {
const int windowSize = 400;
double *analysisWindow = new double[windowSize];
double *synthesisWindow = new double[windowSize];
for (int k = 0; k < windowSize/4; k++) {
analysisWindow[k] = 0;
}
for (int k = windowSize/4; k < windowSize*7/8; k++) {
double x = 2 * M_PI * ((k+0.5)/windowSize - 1.75);
analysisWindow[k] = 2.57392230162633461887-1.58661480271141974718*cos(x)+3.80257516644523141380*sin(x)
-1.93437090055110760822*cos(2*x)-3.27163999159752183488*sin(2*x)+3.26617449847621266201*cos(3*x)
-0.30335261753524439543*sin(3*x)-0.92126091064427817479*cos(4*x)+2.33100177294084742741*sin(4*x)
-1.19953922321306438725*cos(5*x)-1.25098147932225423062*sin(5*x)+0.99132076607048635886*cos(6*x)
-0.34506787787355830410*sin(6*x)-0.04028033685700077582*cos(7*x)+0.55461815542612269425*sin(7*x)
-0.21882110175036428856*cos(8*x)-0.10756484378756643594*sin(8*x)+0.06025986430527170007*cos(9*x)
-0.05777077835678736534*sin(9*x)+0.00920984524892982936*cos(10*x)+0.01501989089735343216*sin(10*x);
}
for (int k = 0; k < windowSize/8; k++) {
analysisWindow[windowSize-1-k] = (1 - analysisWindow[windowSize*3/4-1-k]*analysisWindow[windowSize*3/4+k])/analysisWindow[windowSize/2+k];
}
printf("Analysis window:\n");
for (int k = 0; k < windowSize; k++) {
printf("%d\t%.10f\n", k, analysisWindow[k]);
}
double accu, accu2;
for (int k = 0; k < windowSize; k++) {
accu += k*analysisWindow[k]*analysisWindow[k];
accu2 += analysisWindow[k]*analysisWindow[k];
}
for (int k = 0; k < windowSize; k++) {
synthesisWindow[k] = analysisWindow[windowSize-1-k];
}
printf("\nSynthesis window:\n");
for (int k = 0; k < windowSize; k++) {
printf("%d\t%.10f\n", k, synthesisWindow[k]);
}
printf("Mean of square of analysis window as probability density function:\n%f", accu/accu2);
printf("\nProduct of analysis and synthesis windows:\n");
for (int k = 0; k < windowSize/2; k++) {
printf("%d\t%.10f\n", k, analysisWindow[windowSize/2+k]*synthesisWindow[k]);
}
printf("\nSum of overlapping products of windows:\n");
for (int k = 0; k < windowSize/4; k++) {
printf("%d\t%.10f\n", k, analysisWindow[windowSize/2+k]*synthesisWindow[k]+analysisWindow[windowSize/2+k+windowSize/4]*synthesisWindow[k+windowSize/4]);
}
delete[] analysisWindow;
delete[] synthesisWindow;
}
以及与Kiss FFT和优化库一起使用的优化成本函数的源代码：
class WinProblem : public Opti::Problem {
private:
int numParams;
double *min;
double *max;
kiss_fft_scalar *timeData;
kiss_fft_cpx *freqData;
int smallSize;
int bigSize;
kiss_fftr_cfg smallFFTR;
kiss_fftr_cfg smallIFFTR;
kiss_fftr_cfg bigFFTR;
kiss_fftr_cfg bigIFFTR;
public:
// numParams must be odd
WinProblem(int numParams, int smallSize, int bigSize, double* candidate = NULL) : numParams(numParams), smallSize(smallSize), bigSize(bigSize) {
min = new double[numParams];
max = new double[numParams];
if (candidate != NULL) {
for (int i = 0; i < numParams; i++) {
min[i] = candidate[i]-fabs(candidate[i])*(1.0/65536);
max[i] = candidate[i]+fabs(candidate[i])*(1.0/65536);
}
} else {
for (int i = 0; i < numParams; i++) {
min[i] = -1;
max[i] = 1;
}
}
timeData = new kiss_fft_scalar[bigSize];
freqData = new kiss_fft_cpx[bigSize/2+1];
smallFFTR = kiss_fftr_alloc(smallSize, 0, NULL, NULL);
smallIFFTR = kiss_fftr_alloc(smallSize, 1, NULL, NULL);
bigFFTR = kiss_fftr_alloc(bigSize, 0, NULL, NULL);
bigIFFTR = kiss_fftr_alloc(bigSize, 1, NULL, NULL);
}
double *getMin() {
return min;
}
double *getMax() {
return max;
}
// ___ __ 1
// | \ | | | | | | | / |
// | \ | | | | | | | / |
// | \_ | | | | | | | / |
// | \|__ | | | | | | /| |
// | | -----|_______|___ | | | | / | |
// | | | | ----| | | |/ | |
// --------------------------------x-----------------------x---|---- 0
// 0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 15/16
// |-------------------------------| |-------|
// zeroStarts winStarts
//
// f(x) = 0 if 4/8 < x < 7/8
// f(-x)f(x) + f(-x+1/8)f(x-1/8) = 1 if 0 < x < 1/8
double costFunction(double *params, double compare, int print) {
double penalty = 0;
double accu = params[0]/2;
for (int i = 1; i < numParams; i += 2) {
accu += params[i];
}
if (print) {
printf("%.20f", params[0]/2/accu);
for (int i = 1; i < numParams; i += 2) {
printf("+%.20fcos(%d pi x)", params[i]/accu, (i+1)/2);
printf("+%.20fsin(%d pi x)", params[i+1]/accu, (i+1)/2);
}
printf("\n");
}
if (accu != 0) {
for (int i = 0; i < numParams; i++) {
params[i] /= accu;
}
}
const int zeroStarts = 4; // Normally 4
const int winStarts = 2; // Normally 1
int i = 0;
int j = 0;
freqData[j].r = params[i++];
freqData[j++].i = 0;
for (; i < numParams;) {
freqData[j].r = params[i++];
freqData[j++].i = params[i++];
}
for (; j <= smallSize/2;) {
freqData[j].r = 0;
freqData[j++].i = 0;
}
kiss_fftri(smallIFFTR, freqData, timeData);
double scale = 1.0/timeData[0];
double tilt = 0;
double tilt2 = 0;
for (int i = 2; i < numParams; i += 2) {
if ((i/2)%2) {
tilt2 += (i/2)*params[i]*scale;
} else {
tilt2 -= (i/2)*params[i]*scale;
}
tilt += (i/2)*params[i]*scale;
}
penalty += fabs(tilt);
penalty += fabs(tilt2);
double accu2 = 0;
for (int i = 0; i < smallSize; i++) {
timeData[i] *= scale;
}
penalty += fabs(timeData[zeroStarts*smallSize/8]);
penalty += fabs(timeData[winStarts*smallSize/16]*timeData[smallSize-winStarts*smallSize/16]-0.5);
for (int i = 1; i < winStarts*smallSize/16; i++) {
// Last 16th
timeData[bigSize-winStarts*smallSize/16+i] = timeData[smallSize-winStarts*smallSize/16+i];
accu2 += timeData[bigSize-winStarts*smallSize/16+i]*timeData[bigSize-winStarts*smallSize/16+i];
}
// f(-1/8+i)*f(1/8-i) + f(i)*f(-i) = 1
// => f(-1/8+i) = (1 - f(i)*f(-i))/f(1/8-i)
// => f(-1/16) = (1 - f(1/16)*f(-1/16))/f(1/16)
// = 1/(2 f(1/16))
for (int i = 1; i < winStarts*smallSize/16; i++) {
// 2nd last 16th
timeData[bigSize-winStarts*smallSize/8+i] = (1 - timeData[i]*timeData[bigSize-i])/timeData[winStarts*smallSize/8-i];
accu2 += timeData[bigSize-winStarts*smallSize/8+i]*timeData[bigSize-winStarts*smallSize/8+i];
}
// Between 2nd last and last 16th
timeData[bigSize-winStarts*smallSize/16] = 1/(2*timeData[winStarts*smallSize/16]);
accu2 += timeData[bigSize-winStarts*smallSize/16]*timeData[bigSize-winStarts*smallSize/16];
for (int i = zeroStarts*smallSize/8; i <= bigSize-winStarts*smallSize/8; i++) {
timeData[i] = 0;
}
for (int i = 0; i < zeroStarts*smallSize/8; i++) {
accu2 += timeData[i]*timeData[i];
}
if (print > 1) {
printf("\n");
for (int x = 0; x < bigSize; x++) {
printf("%d,%f\n", x, timeData[x]);
}
}
scale = 1/sqrt(accu2);
if (print) {
printf("sqrt(accu2) = %f\n", sqrt(accu2));
}
double tSpread = 0;
timeData[0] *= scale;
double tMean = 0;
for (int i = 1; i <= zeroStarts*smallSize/8; i++) {
timeData[i] *= scale;
// tSpread += ((double)i)*((double)i)*(timeData[i]*timeData[i]);
double x_0 = timeData[i-1]*timeData[i-1];
double x_1 = timeData[i]*timeData[i];
tSpread += ((double)i)*((double)i)*(x_0 + x_1)*0.5 - ((double)i)*(2.0/3*x_0 + 1.0/3*x_1) + 0.25*x_0 + 1.0/12*x_1;
double slope = timeData[i]-timeData[i-1];
if (slope > 0) {
penalty += slope+1;
}
tMean += x_1*i;
if (timeData[i] < 0) {
penalty -= timeData[i];
}
}
double x_0 = timeData[0]*timeData[0];
for (int i = 1; i <= winStarts*smallSize/8; i++) {
timeData[bigSize-i] *= scale;
double x_1 = timeData[bigSize-i]*timeData[bigSize-i];
tSpread += ((double)i)*((double)i)*(x_0 + x_1)*0.5 - ((double)i)*(2.0/3*x_0 + 1.0/3*x_1) + 0.25*x_0 + 1.0/12*x_1;
x_0 = x_1;
tMean += x_1*(-i);
}
tMean /= smallSize;
penalty += fabs(tMean);
if (tMean > 0) {
penalty += 1;
}
tSpread /= ((double)smallSize)*((double)smallSize);
if (print) {
printf("tSpread = %f\n", tSpread);
}
kiss_fftr(bigFFTR, timeData, freqData);
double fSpread = 0;
x_0 = freqData[0].r*freqData[0].r;
for (int i = 1; i <= bigSize/2; i++) {
double x_1 = freqData[i].r*freqData[i].r+freqData[i].i*freqData[i].i;
fSpread += ((double)i)*((double)i)*(x_0 + x_1)*0.5 - ((double)i)*(2.0/3*x_0 + 1.0/3*x_1) + 0.25*x_0 + 1.0/12*x_1;
x_0 = x_1;
}
if (print > 1) {
for (int i = 0; i <= bigSize/2; i++) {
printf("%d,%f,%f\n", i, freqData[i].r, freqData[i].i);
}
}
fSpread /= bigSize; // Includes kiss_fft scaling
if (print) {
printf("fSpread = %f\n", fSpread);
printf("%f,%f,%f\n", tSpread, fSpread, tSpread*fSpread);
}
return tSpread*fSpread + penalty;
}
double costFunction(double *params, double compare) {
return costFunction(params, compare, false);
}
int getNumDimensions() {
return numParams;
}
~WinProblem() {
delete[] min;
delete[] max;
delete[] timeData;
delete[] freqData;
KISS_FFT_FREE(smallFFTR);
KISS_FFT_FREE(smallIFFTR);
KISS_FFT_FREE(bigFFTR);
KISS_FFT_FREE(bigIFFTR);
}
};
这取决于窗口的上下文。传统上开发的窗口化旨在用于估计功率谱密度的 Blackman-Tukey 方法。这是相关图方法的一般形式，其中利用了离散时间 Wiener-Khinchin 定理。回想一下，这通过离散时间傅里叶变换将自相关序列与功率谱密度联系起来。
因此，窗户的设计考虑了几个标准。首先，它们必须在原点获得统一增益。这是为了保持信号自相关序列中的功率，因为 rxx[0] 可以被认为是样本功率。接下来，窗口应该从原点逐渐变细。这是出于多种原因。首先，为了成为有效的自相关序列，所有其他滞后必须小于或等于原点。其次，这允许较低滞后的较高权重，这是使用大多数样本以极大的信心计算得出的，而较高滞后的权重较小或为零，由于可用于它们的数据样本量的减少而具有增加的方差计算。这最终导致更宽的主瓣和随后降低 PSD 估计的分辨率，
最后，如果窗口具有非负光谱，也是非常需要的。这是因为使用 Blackman-Tukey 方法，您可以将最终估计的偏差视为与窗口频谱卷积的真实功率谱密度。如果此窗口频谱具有负区域，则功率谱密度估计中可能存在负区域。这显然是不希望的，因为它在这种情况下几乎没有物理意义。此外，您会注意到 Blackman-Tukey 方法中没有幅度平方运算。这是因为，通过实数偶数自相关序列乘以实数偶数窗口，离散傅里叶变换也将是实数偶数。在实践中，您会发现通常被量化的非常小的负分量。
由于这些原因，窗口也是奇数长度，因为所有有效的自相关序列也是如此。现在，仍然可以（并且已经完成）的是在周期图方法的上下文中进行窗口化。即对数据加窗，然后取加窗数据的幅度平方。这不等同于 Blackman-Tukey 方法。通过一些统计推导，您可以发现它们的平均行为相似，但一般不相似。例如，在 Welch 或 Bartlett 方法中对每个段使用加窗以减少估计的方差是很常见的。所以本质上，使用这些方法，动机部分相同，但不同。在这些方法中，功率是通过例如划分窗口能量来归一化的，而不是仔细加权窗口滞后。
因此，希望这可以将窗口及其起源以及它们为什么是对称的上下文关联起来。如果您对为什么选择非对称窗口感到好奇，请考虑傅里叶变换的对偶属性的含义，以及卷积功率谱密度估计对您的应用意味着什么。干杯。
加窗的最初目的是确保（由 DFT 假定为周期性）信号在开始时与结束时相比没有尖锐的瞬态。代价是朝向（对称）窗口中心的频率将在随后的 DFT 中得到更多加权和表示。
在所有这些背景下，我可以想象人们会想要使用非对称窗口来强调通过 DFT 分析的信号中的局部时间特征。但是，如果信号的端点在加窗后幅度不大致相同，这可能会以 DFT 期间更宽的波瓣宽度为代价。