Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It only takes a minute to sign up.
Sign up to join this community
Common non-rectangular window functions all seem to be symmetric. Is there ever a case when one would want to use a non-symmetric window function before an FFT? (Say if the data on one side of the FFT aperture were considered a bit more important than data on the other, or less noisy, etc.)
If so, what kinds of asymmetric window functions have been studied, and how would they affect the frequency response compared to a (more lossy?) offset symmetric window?
I'll use the shorthand window for "window function".
With audio, any processing that creates something akin to pre-ringing or pre-echo will sound sloshy like a low bit-rate mp3. This happens when localized energy of a transient or an impulse is spread backwards in time, for example by modification of the spectral data in lapped transforms such as lapped modified discrete cosine transform (MDCT). In such processing, audio is windowed by overlapping analysis windows, transformed, processed in frequency domain (like data-compressed to a smaller bitrate), windowed again with a synthesis window and summed back together. The product of the analysis and synthesis window must be such that the overlapping windows sum to unity.
Traditionally the window functions used have been symmetric, and their width has been a compromise between frequency selectivity (long window) and time-domain artifact avoidance (short window). The wider the window, the more back in time the processing can spread the signal. A more recent solution is to use an asymmetrical window. The two windows used can be mirror images of each other. The analysis window drops from peak to zero fast so that impulses are not "detected" much in advance, and the syntheis window rises from zero to peak fast, so that the effects of any processing are not spread much backwards in time. Another advantage of this is low latency. The asymmetrical windows can have good frequency-selectivity, and can replace variable-sized symmetrical windows in audio compression, like a kind of cure-all. See M. Schnell, M. Schmidt, M. Jander, T. Albert, R. Geiger, V. Ruoppila, P. Ekstrand, M. Lutzky, B. Grill, “MPEG-4 Enhanced Low Delay AAC - a new standard for high quality communication”, 125th AES Convention, San Francisco, CA, USA, preprint 7503, Oct. 2008 and another conference paper where they show also the magnitude of the Fourier transform of their window: Schnell, M., et al. 2007. Enhanced MPEG-4 Low Delay AAC – Low Bitrate High Quality Communication. In 122th AES Convention.
Figure 1. Illustration of use of asymmetrical windows in lapped analysis-processing-synthesis. The product (black dashed) of the analysis-window (blue) and the synthesis window (yellowish orange) sums to unity with the window from the previous frame (grey dashed). Further constraints are needed to guarantee perfect reconstruction when using MDCT.
Discrete Fourier transform (DFT, FFT) could be used instead of MDCT, but will in such contexts give redundant spectral data. Compared to DFT, MDCT gives only half of the spectral data while still enabling perfect reconstruction if suitable windows are chosen.
Here is my own asymmetrical window design (Fig. 2) suitable for lapped analysis-processing-synthesis using DFT but not MDCT with which it does not give perfect reconstruction. The window tries to minimize the product of mean-square time and frequency bandwidths (similarly to the confined Gaussian window) while retaining some potentially useful time-domain properties: nonnegative, unimodal with the peak at "time zero" around which the analysis and synthesis windows are mirror images of each other, function and first derivative continuity, zero-mean when the square of the window function is interpreted as an unnormalized probability density function. The window was optimized using differential evolution.
Figure 2. Left: An asymmetrical analysis window suitable for overlapped analysis-processing-resynthesis together with its time-reversed counterpart synthesis window. Right: Cosine window, with the same latency as the asymmetrical window
Figure 3. Magnitude of the Fourier transforms of the cosine window (blue) and the asymmetrical window (orange) of Fig. 2. The asymmetrical window shows better frequency selectivity.
Here is Octave source code for the plots and for the asymmetrical window. The plotting code comes from Wikimedia Commons. On Linux I recommend installing gnuplot, epstool, pstoedit, transfig first and librsvg2-bin for viewing using 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");
You may want to use only every second sample of the window because it starts and ends at zero. The following C++ code does that for you so you don't get any zero samples except for in one quarter of the window that is zero everywhere. For the analysis window this is the first quarter and for the synthesis window this is the last quarter. The second half of the analysis window should be aligned with the first half of the synthesis window for calculation of their product. The code also tests the mean of the window (as a probability density function), and showcases the flatness of overlapped reconstruction.
#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;
}
And the source code for the optimization cost function to be used with Kiss FFT and an optimization library:
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);
}
};
It depends on the context of windowing. Windowing, as it was traditionally developed, was intended for the Blackman-Tukey method of power spectral density of estimation. This is the general form of the correlogram methods, whereby the discrete-time Wiener-Khinchin theorem is exploited. Recall this relates the autocorrelation sequence to the power spectral density through the discrete time Fourier transform.
Therefore, windows were designed with several criteria in mind. First, they had to have unity gain at the origin. This was to preserve power in the signal's autocorrelation sequence, as rxx[0] can be thought of as the sample power. Next, the window should taper from the origin. This is for a number of reasons. First, in order to be a valid autocorrelation sequence, all other lags must be less than or equal to the origin. Second, this allowed for higher weighting of the lower lags, which have been computed with great confidence using most of the samples, and small or zero weighting of the higher lags, which have increasing variance due to the decreasing amount of data samples available for their calculation. This ultimately results in a wider main lobe and subsequently decreased resolution in the PSD estimate, but a significant reduction in the side lobes that can bias your main lobe of sinusoidal components.
Finally, it's also highly desired if the windows have a nonnegative spectrum. This is because with the Blackman-Tukey method, you can think of the bias of the final estimate as the true power spectral density convolved with the window spectrum. If this window spectrum has negative regions, its possible to have negative regions in your power spectral density estimate. This obviously is undesired, as it has little physical meaning in this context. In addition, you'll note there's no magnitude squared operation in the Blackman-Tukey method. This is because, with a real and even autocorrelation sequence multiplied by a real and even window, the discrete Fourier transform will also be real and even. In practice, you'll find very small negative components that usually are quantized out.
For these reasons, windows are also of odd length because all valid autocorrelation sequences are as well. Now, what can still be done (and is done) is windowing in the context of the periodogram methods. That is, window the data, and then take the magnitude squared of the windowed data. This is not equivalent to the Blackman-Tukey method. You can find, through some statistical derivations, that they behave similarly on average, but not in general. For example, it's quite common to use windowing for each segment in Welch's or Bartlett's method to decrease the variance of the estimates. So in essence, with these methods, the motivation is in part the same, but different. Power is normalized in these methods by dividing out the window energy for example, instead of careful weighting of the window lags.
So, hopefully this contextualizes windows and their origins, and why they are symmetric. If you're curious as to why one may chose an asymmetric window, consider the implications of the duality property of the Fourier transform, and what convolving your power spectral density estimate means for your application. Cheers.
The original point of windowing is to make sure that the (assumed periodic by the DFT) signal has no sharp transients at the beginning compared to the end. The cost is that frequencies towards the center of the (symmetric) window will be more weighted and represented in the subsequent DFT.
With all that in the background, I can imagine that one would want to use an asymmetric window to accentuate local temporal features in the signal being analyzed via the DFT. However this might come at the cost of a wider lobe width during DFT, if the end points of your signal are not roughly the same amplitude after windowing.
