Safe Point Thinning Algorithm - How to? I am relatively new to image processing and am working on a project that requires the thinning of binary images. In the literature, I've seen references to the Safe Point Thinning Algorithm(SPTA) quite a lot, but have failed to find any definitive source to the algorithm, maybe because most of the articles quoted to describe the SPTA are not in the public domain. Can you please explain SPTA for me? Or even point me to (free)resources that do? The pros and cons of the algorithm are welcome. I could also use any alternative algorithm that is "good enough".

SPTA seems to be a quite old skeleton algorithm that is even hard to find behind paywalls.
Here is a link to a skeletonization paper that references it and comes with free access : http://dspace1.acs.ucalgary.ca/bitstream/1880/45551/2/1987-287-35.pdf

But roughly speaking the idea of skeleton algorithms is to look of a simple shape "midway" inside the detected blob. For example, if you have detected a rectangle, then its skeleton is going to be the points inside the rectangle that are at equal distance from its boundaries, which makes a shape like this (the part inside the rectangle):

--------
|>----<|
--------

If it helps, think of the way children draw characters wit just a few strokes...

-- EDIT --

OK, if I had to choose some reference to implement (but I understood from the question that SPTA was actually used, not to be replaced), I would check the papers on this page:
http://www.ensta-paristech.fr/~manzaner/publis.html

Disclaimer: I'm not this guy, I didn't write the papers, but I would start here because:

the papers are available for download
A. Manzanera worked a lot on limited hardware (e.g. embedded stuff), so I guess the algorithms are efficient and not so tricky to implement


DTFT of inverse of any function In my book solution is given like this.
But i am solving like this , am i doing wrong??

Yes, you're doing it wrong. :) Inverse systems in time domain satisfy

$$h[n] * h_i[n] = \delta[n]$$

where $*$ denotes convolution and $h_i[n]$ is the impulse response of the inverse system. If you DTFT this equation, you'll get a).


Direct path and reflections of ultrasound transmission I am using the following system:
Ultrasonic chirp signals with central frequency 100kHz are being transmitted using ultrasonic sensor. After the transmission through the medium(air/solid/liquid/combination) they are being received with hydrophone. The system is being controled by DSP and signal generation and acquisition is obtained via DAC and ADC. Our aim is to detect the time of arrival between transmission and reception.
My mentor suggested that I could use something with FFT to observe reflections, but I am struggling to find information on the internet about it. Does anyone have a suggestion how that could be done using FFT, or any other method?

That sounds like a usual FMCW sonar.
It's technically identical to FMCW radar. So, multiply your receive with your transmit signal, you get a signal at the beat frequency, which is proportional to the range.
Standard FMCW approaches; there's a wealth of info out there.


Ambiguity in the IFFT process in OFDM I am still trying to iron out some ambiguities in my understanding of the IFFT component of OFDM modulation schemes.

So here we have a QAM symbol $s_0$ being multiplied with the subcarrier for that frequency bin, $s_1$ being multiplied with carrier 1, $s_2$ carrier 2, and it shows clearly that the sum of these products produces the OFDM symbol in the time domain. This would logically follow the equation:

$
c\left(t\right) = \sum^{N-1}_{n=0}{x_n e^{-i\frac{2\pi n}{N}}}
$

Where c(t) is the output symbol, x is the list of complex QAM symbols and n is the subcarrier index, N is the number of subcarriers. (Note there is no 1/N at the start and i'm not sure why that needs to be introduced)

Why is it that IDFT is stated on Wikipedia to be

$
x_n = \frac{1}{N} \sum^{N-1}_{k=0}{X_k . e^{i 2 \pi k n / N}}
$

Where x is the output list and xn is merely a single output sampling point in the time domain.

Here is an example of application of the DFT/IDFT formula :

Let :
$N = 4$
and :
$
x = \left(\begin{array}{c}x_0\\ x_1\\ x_2\\ x_3\\
\end{array}
\right)
=
\left(\begin{array}{c}1\\ 2-i\\ -i\\ -1+2i\\
\end{array}
\right)
$

then

$X_0 = e^{-i 2 \pi 0.0/4} + e^{-i 2 \pi 0.1/4}.\left(2-i\right) + e^{-i 2 \pi 0.2/4}.\left(-i\right)+ e^{-i 2 \pi 0.3/4}.\left(-1+2i\right) = 2$
$X_1 = e^{-i 2 \pi 1.0/4} + e^{-i 2 \pi 1.1/4}.\left(2-i\right) + e^{-i 2 \pi 1.2/4}.\left(-i\right)+ e^{-i 2 \pi 1.3/4}.\left(-1+2i\right) = -2-2i$
$X_2 = e^{-i 2 \pi 2.0/4} + e^{-i 2 \pi 2.1/4}.\left(2-i\right) + e^{-i 2 \pi 2.2/4}.\left(-i\right)+ e^{-i 2 \pi 2.3/4}.\left(-1+2i\right) = -2i$
$X_3 = e^{-i 2 \pi 3.0/4} + e^{-i 2 \pi 3.1/4}.\left(2-i\right) + e^{-i 2 \pi 3.2/4}.\left(-i\right)+ e^{-i 2 \pi 3.3/4}.\left(-1+2i\right) = 4+4i$

$
X =
\left(\begin{array}{c}X_0\\ X_1\\ X_2\\ X_3\\
\end{array}
\right)
=
\left(\begin{array}{c}2\\ -2-2i\\ -2i\\ 4+4i\\
\end{array}
\right)
$

If we pretend this worked example is IDFT and not DFT, it clearly shows that for each output sampling point, Xk, in the time domain, every point in the constellation was multiplied by a single frequency bin and then summated. This is directly backed up by the fact that IDFT is stated to have a complexity of O(n^2) (as far as I'm aware IFDT and IFFT are the same equation but implemented differently algorithmically).

Why is this? This directly contradicts many OFDM examples and diagrams I've seen and it fundamentally makes sense to me that a single frequency would be phase and amplitude shifted by a single complex input from the constellation, and not by all of them.

The answers to this are now on my site, which I'm not going to keep updating into this answer when I modify it


Linear predictive model convolution Accidentally asked this question in the general area and was told to ask here, so...

I've been trying to develop a lightweight, relatively-fast-to-decode sound compression format for use in my gaming projects (perfect reproduction is not needed so I only use 16 bit data).

The idea is to split the sound data into 14-sample frames and use linear prediction to reduce the data size by only storing the residuals. To make it even lighter, the residuals are then quantized to 4 bits per sample by reducing their precision so that they are stored as scaled residuals, where the scale is dictated by the frame header. To make the signal less noisy, 16 linear prediction models are generated that best suit the signal in the file.

Each frame ends up being 64 bits (4 bits for the LP model selection, 4 bits for the residual scale (2^n, so only n is stored), and 14*4 bits for the residual data).

What I've been doing is simply solving the Yule-Walker equations for each frame to get the coefficients and then saving them in an array. After calculating the coefficients for all the frames, I then try to create 16 coefficients based on the Euclidean distance from each LP model and averaging them out. That is to say...

Zero-out target buffer

For i = 1, i <= 16, i++
  Create and zero-out a temporary buffer for i LP models

  For each saved LP model
    Find the model in the original target buffer [up to i] with the smallest Euclidean distance between it and this model
    Accumulate the coefficients to the temporary buffer

  Average out the temporary buffer and store it to the final LP model array [again up to i]

As expected, this doesn't yield very good results at all if only because I imagine it's similar to uniform colour quantization - it has very little relation to the actual data.

After Googling for a few months, the only thing I can come up with is 'sparse convolution' but I'm not entirely sure what is meant by that when it comes to LPC or if it's even what I'm after.

Given these parameters (that is: generate 16 LPC models that will minimize the residual error given that the frames are 14 samples long), how would you go about doing it?

Averaging AR-coefficients is a risky business - this could cause drastic shifts in the position of the zeroes of the corresponding polynomial, resulting in instabilities!

If I were to approach the problem the first thing I would try would be the following:

Compute the spectrogram (short-term fourier transform)
Run VQ on the spectral slices with a codebook size of 16, using euclidian distance on the log
Once you have 16 spectral template, convert them back to autocorrelation functions, solve yule-walker for the 16 templates and you get your 16 models.

Or another "safe" way:

For each frame, extract the AR coefficients.
Convert them into reflection coefficients or any other alternative representation like line spectral pairs which can be safely interpolated.
Run VQ on these with a codebook size of 16 to get your 16 models.
Convert back to AR coefficients.

But I really doubt this is a path worth exploring. My two intuitions would be that:

If an audio file is long and diverse enough the 16 specialized models you would get out of it won't be very different form the 16 "universal" models you would obtain out of running the procedure once and for all on a large collection of audio files.
Learning 16 models from a large collection of audio files is likely to give you the kind of naive "fixed predictors" used in FLAC (polynomial interpolators).

One related result is that running a KLT on mel-frequency spectra on a sufficiently long audio clip is going to yield something very similar to the DCT base - that's why audio codecs like mp3 or AAC are not using per-file optimized transforms - the DCT is a good, universal one, and the advantage of using universal transforms is that they don't even require any side-information. This is what motivates my intuition - you are looking at a large dataset (an entire audio file), 16 is small. You'd rather think of 16 good universal predictors and use them than instead of trying to optimize your 16 predictors per file.


Why is the first value in digital first-order IIR filter impulse response not the largest? When I create a digital first-order IIR low-pass filter with scipy (code below), I get the following coefficients:

B: [ 0.1367  0.1367]
A: [ 1.     -0.7265]

The impulse response of this filter does not look like most text-book examples in that the first value is not the largest.

Why is this so?  Does it have to be so?

If I want to make an FIR filter to mimic such a first-order IIR low-pass filter, should it be a normal exponential decay shape, or does it need the first value "tweaked" in this manner?

import numpy as np
from scipy.signal import butter, lfilter
import matplotlib.pyplot as plt

b, a  = butter(1, .1, btype='lowpass')
print('B:', b)
print('A:', a)

x = np.zeros(1000)
x[0] = 1.

xf = lfilter(b, a, x)

plt.figure()
plt.plot(x,'o', label='impulse')
plt.plot(xf,'s-', label='response')
plt.xlim(-1, 10)
plt.ylim(-.1, .3)
plt.legend()

A causal first-order IIR filter is characterized by the following difference equation:

$$y[n]=b_0x[n]+b_1x[n-1]-a_1y[n-1]\tag{1}$$

with $x[n]$ the input signal, and $y[n]$ the output signal. The impulse response of that system can be computed via the $\mathcal{Z}$-transform or otherwise, and it turns out to be

$$h[n]=b_0\delta[n]+(-a_1)^{n-1}(b_1-b_0a_1)u[n-1]\tag{2}$$

where $\delta[n]$ is the unit impulse, and $u[n]$ is the unit step sequence. Assuming that the filter is stable, we can conclude that $|a_1|<1$ is satisfied. Consequently, $|h[0]|>|h[n]|$ for $n>0$ is satisfied if and only if

$$|b_0|>|b_1-b_0a_1|\tag{3}$$

holds.

With $b_0=b_1$ and $a_1<0$ as in your example, condition $(3)$ can never be satisfied, and, consequently, the maximum magnitude of the impulse response is not assumed at index $n=0$.

EDIT:

I add more details referring to your filter design problem. A Butterworth low pass filter's frequency response satisfies $H(1)=1$ (DC) and $H(-1)=0$ (Nyquist). With

$$H(z)=\frac{b_0+b_1z^{-1}}{1+a_1z^{-1}}\tag{4}$$

we have the following requirements on the coefficients:

$$H(1)=\frac{b_0+b_1}{1+a_1}=1\\
H(-1)=\frac{b_0-b_1}{1-a_1}=0\tag{5}$$

which gives

$$b_0=b_1\quad\textrm{ and }\quad a_1=2b_0-1\tag{6}$$

From the stability condition $|a_1|<1$ we get

$$|2b_0-1|<1\Longrightarrow 0<b_0<1\tag{7}$$

Combining the constraints $(6)$ and $(7)$ with the inequality $(3)$ yields

$$1 >b_0>\frac12\tag{8}$$

If $(8)$ is satisfied, $h[0]$ will have the largest magnitude, otherwise it won't. Now it turns out that the exact value of $b_0$ depends on the desired cut-off frequency of the low pass filter. For cut-off frequencies less than $\frac12$ (i.e., less than half the Nyquist frequency), you get $b_0<\frac12$, and hence $|h[0]|$ is smaller than the maximum magnitude of the impulse response. Only for cut-off frequencies greater than half the Nyquist frequency do you get $|h[0]|>|h[n]|$, $n>0$.

So to directly answer the question in your title: the first impulse response value is not the largest (in magnitude) because the filter coefficients do not satisfy inequality $(3)$, and in your example of a (first-order) Butterworth low pass filter this happens because the cut-off frequency is smaller than half the Nyquist frequency.


FFT with asymmetric windowing? 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);
  }
};


Sending a Filtered vs Non-Filtered speech over a Communication Link To send a $10 \, KHz$ bandwidth of speech (Non Filtered speech), Nyquist criteria says that sampling duration should be $50 \, \mu s$. While for a $4 \, KHz$ bandwidth of filtered speech, sampling duration should be $125 \, \mu s$ (Much larger than $50 \, \mu s$.).

Now an excerpt taken from my textbook:

"In general, the wider the $BW$ of a signal, the more "expense" it will be to transmit the signal. So it is much "cheaper" to send a filtered speech over a communication link."

What are the factors related to the words "expense" and "cheaper" mentioned here? What does author mean by using these two words?

It can mean a number of things; some examples are:

Converters (analog to digital as well as digital to analog) go up in price and in power consumption as the sampling rate increases.
In very general terms, analog circuits capable of handling signals of larger bandwidth are more expensive than narrower-band circuits. Compare the cost of a 741 op-amp to a 2082's.
In the digital domain, a higher sampling rate implies more calculations per second, which in turn implies higher-performance digital circuits and more power consumption.
Channel bandwidth has a cost. In wireless communications, you either have to license a band or use an unlicensed band (such as ISM). ISM bands are subject to interference and put limits on the power you can radiate. In wireline communications, wire cost goes up with bandwidth -- compare the cost of twisted pair to coax to thunderbolt wires.


Does MatLab's `ftt` function gives a Fourier Transformed signal in terms of frequency $f$ or angular frequency $\omega$? When we use MatLab's fft function to Fourier transform a time signal, do we obtain a Fourier space signal in terms of frequency $f$, or angular frequency $\omega$?

Let say we have

tDelta = 0.1;
tMax = 8;
t = 0:tDelta:tMax;

sigFreq = 3*(2*pi);
s = sin(sigFreq*t);
s_FT = fft(s);

Here's what I think is going on, please correct me if I'm wrong:

It seems that tDelta = 0.1 seconds means that the sampling rate is 1/tDelta = 10 Hz.
And then when we Fourier transform the signal $s(t)$ with ftt MatLab gives us the spectrum $\hat{s}(f)$ with the frequency $f\in[-5,5]$ Hz.
To convert this to angular frequency $\omega$ we must multiple $f$ by $2\pi$.
So finally, in terms of angular frequency we have $\hat{s}(\omega)$ with $\omega\in[-5*(2\pi),5*(2\pi)]$.

Is this correct or have I have misinterpreted how fft works? I'm really not certain about when we have are dealing with frequency $f$ and when we are dealing with angular frequency $f$?

There are two things to clarify here before the answer becomes self-evident.

The trigonometric functions $\cos, \sin, \dots$
The relationship between phase $\phi$ and frequency $f$.

The trigonometric functions $\cos(\phi)$ and $\sin(\phi)$ trace the periphery of a circle. Where we are on the circle is denoted by the argument of the function \phi. But, in general, once we are around the circle once, the values keep repeating periodically.

Now, suppose that we say $\cos(4)$ ... Where is that on the circle? Is this on the "4" or is it further down its "multiples" sequence?

If we said, $\cos(2 \pi 4)$, then we know this represents 4 times around a circle of $2 \pi$. If we say $\cos(2 \pi 0.5)$ then we know that's halfway around a full circle.

Please note that so far we are talking about instantaneous values, akin to saying "What is the cos of 60 degrees? What is the cos of (360 + 60) degrees? What is the cos of (2*360 + 60) degrees? What is the cos of (k*360 + 60) degrees?" and so on.

Let's bring time in to this. Let's now say, not "Where am I, around the circle?" but "How fast am I going around the circle?".

$$ y = \cos(2 \pi t) $$

Where $t$ is time and time is a variable we cannot stop, it starts (usually) at $0.0$ at some arbitrary set off point (i.e., $t=0$ when we power on a system) and increases for ever beyond our control. It also has units of seconds.

In $1.0$ second, we have been around the $2 \pi$ circle once, in $2.0$ seconds we have been around it twice...and so on.

This represents periodic events with a period of $1.0$ sec, the time it takes to go around the circle.

What if we have things that happen more frequently than once every second?

This is easy to deduct, we just have to go around the circle faster. How much faster? I don't care, some factor $f$ faster than 1 rev per second.

$$ y = cos(2 \pi t f + \theta )$$

So, now, the phase at which we find ourselves moving around the circle is $\phi = 2 \pi t f$ plus some $\theta$ term to account for any starting point around the circle.

If we differentiate phase with respect to time we get frequency. That is, the rate at which we move around the circle.

The rate has a $2 \pi$ factor, yes, but that offsets everything and means number of times around the circle. The key factor here is $f$.

So, $2 \pi 440$ radians per second is another way of saying 440 Hz.

Let's now turn to the Fast Fourier Transform (FFT).

We have a system that samples at an $Fs = 44.1$ Hz. We use it to acquire a buffer of 441 samples which we then pass through the FFT. The FFT returns a result that contains 441 complex values, symmetrical about the middle point. The question is, what is the frequency of "bin" 44?

It is $44 \times \frac{44100}{441} = 4400$ Hz. And in this formula we used Hz throughout.

Hope this helps.


What does the word Entropy refer to in Entropy encoding I have come know of Huffman and Arithmatic encoding that try to curtail redundant data and thus compress data. Why are they called Entropy encoding? What does that mean? Why not just call them compression encoding?

The term entropy comes from the field of information theory.  Information theory posits that unexpected messages bear more information than expected ones.  For instance, saying that there will be a solar eclipse tomorrow has much more information (assuming it's true) than saying the sun will rise tomorrow.  Saying the sun will rise tomorrow, while true (and important!), has very little information for most recipients because they have pretty much 100% expectation of it.

Anyway, you can prove mathematically that information is maximized when all of the possible messages have equal probability.  This maximizes the randomness or "entropy".

A perfect, lossless compression scheme could theoretically reduce the number of bits of a message to $ceil(log_2(entropy))$ but no further.  Lossy compression schemes can reduce further, but information will be lost.

I personally am not fond of the term "entropy" because I associate it more with the 2nd law of thermodynamics.  I believe that that association was intentional though, in the sense that entropy (in the thermodynamics sense) implies randomness and chaos, while low entropy implies orderliness.


How do we compute distrubtions of the value of a random process conditional on initial conditions? Suppose I have a stationary process $\phi(t)$ with a known autocorrelation function

$$ A(\tau) \equiv \langle \phi(0) \phi(\tau) \rangle$$

and suppose I also know that $\phi(t)$ is Gaussian distributed.

If I have a particular realization of the process in which $\phi(0)=\phi_0$, what is the conditional distribution of $\phi(t)$ for later times in that realization?

"I know that $\phi(t)$ is Gaussian distributed" is not the same as saying
"$\{\phi(t)\}$ is a Gaussian process" but I will assume that the latter is meant. With this assumption, and the additional assumption that the process is
wide-sense-stationary (the autocorrelation function is listed as having only one argument), the process is strictly stationary and thus
$\phi(t)$ and $\phi(0)$ are jointly Gaussian random variables. They have
the same mean $\mu$ and variance $\sigma^2$, and
their correlation coefficient is $\rho$ where
$$\mu = \sqrt{\lim_{t\to \infty} A(t)}, \qquad \sigma^2 = A(0)-\mu^2,
\qquad \rho(t) = \frac{A(t)-\mu^2}{\sigma^2}.$$
Consequently the conditional distribution of $\phi(t)$ given
that $\phi(0)$ has taken on value $\phi_0$ is Gaussian
with mean $\mu + \rho(t)(\phi_0 - \mu) = \rho(t) \phi_0 + (1-\rho(t))\mu$
and variance $\sigma^2(1-\rho(t)^2)$. Note that as $t\to\infty$,
the conditional distribution of $\phi(t)$ approaches the
unconditional distribution of $\phi(t)$: the distant past affects
the present state less and less.

For more information about Gaussian processes, see, for example,
this answer to a
different question.


Approximation of Hilbert Transform Using Very Short Hilbert Transform FIR / IIR Filter Hilbert transform is a quite sensitive topic here, since Gabor's paper,  Theory of Communication,  J. Inst. Electr. Engineering, London, 1946. Perhaps even more important than the Fourier transform.

There are short Gaussian FIR approximations, smoothing filters, or polynomial fits.

What would be methods for very short (say, less than $16$ taps) Hilbert tranform approximations of a signal, possibly including a weighting factor? I am intending to use them on small chunks of signals, in moving frames that progress along data as new samples are acquired. Signals can be considered freed from a potential low-frequency trend (e.g. polynomial).

Fey features would be:

a management or reduction of small-sized frame artifacts (I am aware of those for short Fourier frames, much less with short Hilbert transforms),
potential asymmetric windowing to add a forgeting factor to older samples,
a recursive formulation, in the spirit of the sliding DFT, would be a plus.

Finally, an IIR formulation, such as the one of the recursive least squares (RLS) adaptive filters, or the exponentially-weighted moving average filter, would be very interesting as well.

This is more like an extended comment to chart the possible answers.

Hilbert transform is a frequency domain 90-degree phase shift of the signal. It has an antisymmetrical impulse response around time = 0. You specify that the approximation shall be causal (EDIT: this requirement has since been removed), so I think you need to reference the phase shift to that provided by another filter. If we denote the two filters F1 and F2, some usable alternatives are:

Case | Filter F1    | Filter F2           | Phase diff     | Magnitude frequency responses
-----+--------------+---------------------+----------------+------------------------------
  C1 | pure delay   | antisymmetrical FIR | exactly 90 deg | F2 has ripple
  C2 | pure delay   | all-pass IIR        | approx 90 deg  | No ripple
  C3 | all-pass IIR | all-pass IIR        | approx 90 deg  | No ripple

In all those cases you have two filters whose phase frequency responses have a phase difference of exactly or close to 90 degrees over the desired frequency band.

If the filter outputs are used "in quadrature" (summed, with the other multiplied by the imaginary unit) to create an analytic signal, any phase frequency response ripple will become magnitude frequency response ripple of the composite negative frequency removal filter.


Generating a sound where end frequency is 1/2 of the start frequency I am very very new to signals, and I am trying to generate sounds based on some parameters from a database. So, the database stores the following:

startFrequency
endFrequency
startAmplitude
endAmplitude
duration

I am generating a WAV file based on this information using the following code:

...
double currentFrequency;
double currentAmplitude;

int samplerate = 44100;
int samples = ceil(duration * samplerate);
int* mySampleBuffer = new int[samples];

for(int i=0; i<samples; i++)
{
    currentFrequency = Interpolate(i, 0, samples, startFrequency, endFrequency);
    currentAmplitude = Interpolate(i, 0, samples, startAmplitude, endAmplitude);
    mySampleBuffer[i] = static_cast<int>(INT_MAX * (sin(((double)i * currentFrequency * PI * 2)/(double)samplerate) * (currentAmplitude/maxAmplitude)));
}

writeWAVData("mySound.wav", mySampleBuffer, mySampleBufferSize, samplerate, 1);
delete[] mySampleBuffer;
...

When the currentFrequency reaches 1/2 of the startFrequency, I notice that the sound wave flips (not sure what the correct term for this is). It doesn't seem like this is correct. Could someone help me understand what is going on.

Here is a screenshot of Audacity

For this image, startFrequency is 3380 Hz, endFrequency = 1380 Hz, duration is 0.512 s, and the amplitude is constant from start to finish. I have calculated that at 0.43264 s, the currentFrequency is exactly 1/2 of the startFrequency. In the picture you can see the "flip" happen at about 0.43 seconds.

Think of frequency in this context more in terms of angular velocity. If that velocity is constant, nothing changes from the usual picture. However, if the frequency changes, it is more correct to compute the actual angle that is input into the trig functions as an anti-derivative of the frequency function.

Here the interpolated frequency is

$$f(t)=(1-\tfrac tT)\,f_0+\tfrac tT\,f_T,$$

so the angle at time $t$ is

$$\phi(t)=2\pi\,\int_0^t f(s)\,ds=2\pi\, [t(1-\tfrac{t}{2T})\,f_0+\tfrac{t^2}{2T}]=\pi t\,[(2-\tfrac tT)f_0+\tfrac tT\,f_T]$$

In contrast, your formula is

$$\phi_o(t)=2\pi\,t\,f(t)=2\pi t\,[(1-\tfrac tT)f_0+\tfrac tT\,f_T]$$

has a derivative or local angular velocity of

$$\phi_o'(t)=2\pi\,[(1-2\tfrac tT)f_0+2\tfrac tT f_T]$$

which has a zero at $t=\frac T2\,\frac{f_0}{f_0-f_T}$, which given your data is located at $t=0.43264s$.

Note that the correct formula also reaches the local angular frequency of zero, but at the double of that time interval, so outside $[0,T]$.

Note on implementation: The last form of $\phi(t)$ can be constructed as

 pi*(i*dt)*(f0+interpolate(i,f0,fT,N))

with N=samples and dt=T/N, provided that

interpolate(0,f0,fT,N)==f0 and interpolate(N,f0,fT,N)==fT


Filter wrap-around to avoid dealing with settling time? Using a Butterworth I've been using a Butterworth high pass filter to correct the low frequency distortion of an acceleration waveform integrated to a velocity waveform.  I throw out a small number of the first data points of the velocity waveform due to the filter settling time.

Is it acceptable or even possible to wrap the filter back around to the beginning and process those 'settling' points as well so I don't have to throw out any data?

Maybe you can reverse the direction of the filter when you hit your known settling point and process back to the beginning then reverse again and go all the way to the end?

I'm not sure if this is an obvious or appropriate question, I have no background in DSP.

If your input waveform is periodic, then wrap-around is appropriate. The most straightforward way to do this is: Initialize the filter state to correspond to an all-zeros input history. Start filtering from the beginning of the input data, writing the output to an output buffer (or in-place to the same buffer, but then the filter has to store its state internally). When the time index runs out of the data buffer, wrap it back to the beginning and keep inputting the filter with zeros for some time while adding the filter's output to successive samples of the output buffer.

You can continue this until the absolute values of the filter's state variables all go under some small threshold value decided by you. This gives an upper limit to residual error in the output. The state of a typical infinite impulse response (IIR) filter is equivalent to (or is directly) its n most recent output values, with n depending on the filter's order. You can keep count of how many of the most recent absolute output values were below the threshold, and when the count reaches n you are done. Or simply decide beforehand how long you will continue, as, for a particular filter, that will determine an upper limit to the residual error proportional to the input.

Note: the period of the input waveform should be an integer multiple of the sampling period.


Squared magnitude of the Z-transform I am basically new to the $z$-transform and there are some points regarding its square magnitude that I do not understand.
Basically I do not understand how in slide 4 of PDF, they arrive at the following expression of $|H(z)|^2$:
$$P(z)=|H(z)|^2=H(z)\ast H^*(1/z^*)$$
If we start with
$$p(n)=h(n)\ast h(-n)$$
but directly take the $z$-transform of the above expression (without substituting $h(n)$ by $h^*(-n)$ as they do in PDF) we arrive at:
$$P(z)=H(z)\ast H(1/z)$$
where I am using the following properties of the $z$-transform:
$$x[n] \longleftrightarrow X(z)$$
$$x[-n] \longleftrightarrow X(1/z)$$
However, it is has been said before (StackExchange) that this expression ($P(z)=H(z)\ast H(1/z)$) is only valid for the unit circle, however I fail to see in which part of the derivation I made the imposition that $|z|=1$.
For reference I attach a capture of the lines I am referring to in PDF:

I understand your confusion. What we are looking for is a function of $z$ which equals the squared magnitude $|H(e^{j\omega})|^2$ of a Fourier transform $H(e^{j\omega})=\mathcal{F}\{h[n]\}$ when evaluated on the unit circle $|z|=1$. Note that the inverse discrete-time Fourier transform of $|H(e^{j\omega})|^2$ is given by the convolution $h[n]\star h^*[-n]$, where $^*$ denotes complex conjugation. Since the $\mathcal{Z}$-transform of $h^*[-n]$ is given by $H^*(1/z^*)$ we see that the $\mathcal{Z}$-transform of $h[n]\star h^*[-n]$ is given by
$$\mathcal{Z}\big\{h[n]\star h^*[-n]\big\}=H(z)H^*\left(\frac{1}{z^*}\right)\tag{1}$$
If $h[n]$ is real-valued we obviously have $h^*[-n]=h[-n]$, and $(1)$ is equivalent to
$$\mathcal{Z}\big\{h[n]\star h[-n]\big\}=H(z)H\left(\frac{1}{z}\right)\tag{2}$$
It's important to understand that the equations
$$\big|H(z)\big|^2=H(z)H^*\left(\frac{1}{z^*}\right)\tag{3}$$
and, for real-valued $h[n]$,
$$\big|H(z)\big|^2=H(z)H\left(\frac{1}{z}\right)\tag{4}$$
are only valid on the unit circle $|z|=1$.
In general we of course have
$$\big|H(z)\big|^2=H(z)H^*(z)\tag{5}$$
but $(5)$ is no valid $\mathcal{Z}$-transform of any time domain sequence, because it is not analytic. On the other hand, $(3)$ and $(4)$ are valid $\mathcal{Z}$-transforms, namely of the sequence $h[n]\star h^*[-n]$. On the unit circle they equal the squared magnitude of the Fourier transform of $h[n]$, i.e.,
$$H(z)H^*\left(\frac{1}{z^*}\right)=\big|H(e^{j\omega})\big|^2,\qquad |z|=1\tag{6}$$
Let's look at a simple example. With $h[n]=\delta[n]+\delta[n-1]$, we have $H(z)=1+z^{-1}$. Since $h[n]$ is real-valued we can use Eq. $(4)$:
$$H(z)H\left(\frac{1}{z}\right)=(1+z^{-1})(1+z)=z+2+z^{-1}\tag{7}$$
The actual squared magnitude of $H(z)$ is
$$\big|H(z)\big|^2=|1+z^{-1}|^2=1+2\textrm{Re}\{z^{-1}\}+|z|^{-2}\tag{8}$$
On the unit circle $z=e^{j\omega}$ both $(7)$ and $(8)$ equal
$$\big|H(e^{j\omega})\big|^2=2+2\cos(\omega)\tag{9}$$
However, the inverse $\mathcal{Z}$-transform of $(7)$ equals the inverse Fourier transform of $|H(e^{j\omega})|^2$, whereas $(8)$ is no $\mathcal{Z}$-transform.


Matching headphone responses First post, and this isn't my field of study, but I'm giving this a go anyway:

My latest endeavor is to use a bit of DSP magic to match up various pairs of headphones, at least in terms of their impulse response measured at the entrance to the ear canal. So far I've taken 64 measurements on 3 sets of headphones, in this way:

1) Put on binaural mics

2) Put on headphones

3) Take IR measurements

These are in a nested loop with 4 iterations each, so I end up with 64 measurements. The IRs I get using a program called Aliki, free software in Linux that uses sine sweeps.

My goal is to generate two filters for each headphone pairing (1<->2, 1<->3, 2<->3). Currently I'm eyeballing a "best" IR for each headphone using spectrograms, matching up the position of peaks, padding, then doing FFT division to get a transfer function (which I have read can be a really bad way of doing things). Results are actually not bad sounding, but I wonder if I can do better with a bit more legwork on the data.

So my question is: What's the best way to combine these IRs to get (possibly) better filters?

Thank you.

Why not to simply match the frequency responses and sensitivities? That's what labs are actually doing.

Since you have the IR's for each device, simply calculate the Frequency response of each repetition, average and compare between devices. In Python that is very simple task.


Why use $\chi^2$ test to determine the presence of white noise? I want to test for the presence of broadband noise in a snapshot 1000 complex baseband samples recorded by a software defined radio.

As a follow-up to this post, why was the $\chi^2$ test used? How many degrees of freedom should be used?

Also, how would one extend this approach to complex baseband data? I would assume that I and Q are iid Gaussian random variables. The magnitude of the complex data would then be Rayleigh distributed not Gaussian. Is there a generalization of the $\chi^2$ test for Rayleigh random variables? Or, would I just pick I or Q to operate on?

Update: I was able to find a paper: A test for whiteness. The author outlines a similar process.

If I and Q are normal, then their magnitude square is $\chi^2$ distributed.

Hence, if you check the distribution of the instantaneous power against that distribution, you get information about how much the complex baseband signal is circularly complex normal.

The autocorrelation function has power of the signal as its value for zero shift. For whiteness, the sample autocorrelation approaches an accordingly weighted Dirac impulse with growing observation length.


Relationship between Discrete Deconvolution and Toeplitz Matrices I have 2 vectors, $a$ & $c$, both of length M. I know they are related by $a*b=c$. My goal is to recover $b$. Obviously $b=$deconv$(c,a)$. I am only interested in the first M elements of the convolution (ie discarding the last M-1 elements). I have two methods of doing this deconvolution in MATLAB and I get different answers. I want to know why:

Method1: using MATLABs default algorithm: b= deconv([c;zeros(M-1,1)],a);

Method2: Using toeplitz matrices:

  A=toeplitz(a,[a(1);zeros(M-1,1)]);
  b=A^-1*c

where A is the toeplitz matrix form of the vector a with 0s on upper triangle.

I believe in theory these methods should give me identical answers, but the answers are somewhat different (but still close). Why is this?

Thanks!

I think you are constructing the convolution matrix incorrectly. Let's do an easy example with some numbers to see how it is constructed

Let's imagine that youhave the following vectors of dimensions $M$ and $N$, And let's for the sake of it, say that $M = 5$ and $N = 4$

$\mathbf{a} = [a_{1},a_{2},a_{3},a_{4}]^{T}$

$\mathbf{b} = [b_{1},b_{2},b_{3},b_{4},b_{5}]^{T}$

$
\mathbf{c} = \mathbf{a}*\mathbf{b}
$

In matrix form, if we assume $\mathbf{a}$ is the system response, the convolution is given by

\begin{equation}
    c = \mathbf{A}\mathbf{b}
\end{equation}

with

\begin{equation}
    \mathbf{A} = \left[
    \begin{array}{ccccc}
        a_{1} & 0 & 0 & 0 & 0 \\
        a_{2} & a_{1} & 0 &  0 & 0 \\
        a_{3} & a_{2} & a_{1} &  0 & 0 \\
        a_{4} & a_{3} & a_{2} & a_{1} &  0 \\
        0 & a_{4} & a_{3} & a_{2} & a_{1} \\
        0 & 0 & a_{4} & a_{3} & a_{2} \\
        0 & 0 & 0 & a_{4} & a_{3} \\
        0 & 0 & 0 & 0 & a_{4} \\
    \end{array}\right]
\end{equation}

By doing the convolution with this matrix, you get the full operation including the transients. So the correct matlab code would be

A = toeplitz([a,zeros(1,M-1)],[a(1),zeros(1,M-1)])

(Check help toepliz to verify how it constructs the matrix)

Use this matrix and the result of the convolution is exactly as with the matrix multiplication.

Now, a word of wisdom, as you can clearly see, the full convolution matrix with the transients is not square, so in order to "invert" it and do deconvolution you will need to use the Moore-Penrose Pseudoinverse. And you will get some differences in the beginning and end of the resulting vector, but it is the best you can do, as using the pseudoinverse would be the same as performing least squares, since the deconvolution is given by

$\mathbf{b} = \left(\mathbf{A}^{T}\mathbf{A}\right)^{-1}\mathbf{A}^{T}\mathbf{c}$

This result will always hold, because the dimensions of $A$ are always $(M+N-1)\times M$, so the matrix is always "skinny" and the equation holds.

Try this, it should give you the result you were looking for.


OCR (optical character recognition) for Indian languages - segmentation Can any one help me in how to extract a table from the whole document image and segment it? I have done segmentation for the document image which has only text by using the projection profile. So please help me how to extract the table from the whole document image. I want only table to be extracted from the document. Please give me an idea how to do it. I want to continue my project using the projection profiles only.

i have done the segmentation manually for the document image which was consisting only text using projection profile manually i din't use the tables in it. now i want to extract the contents of the table in an image without using the formulas or any transforms can you help me.

Why do not use hough transform for finding lines and then finding table region?

you can use hough transform to find horizontal and vertical lines. and then extract region of lines.


What's the meaning of the continuity in spectrum analysis? We all know that for any, suitable, kind of signal $f(t)$,  there corresponds a Fourier transform function $F(j\omega)$ such that
$$
f(t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} F(j\omega) e^{i\omega t} d\omega
$$

So, $F(j\omega)$ can be recognized as the component of $f(t)$ along $e^{i \omega t}$. This is similar to orthogonal analysis. Like a cartesian vector $\vec{a} = (1,2,3)$ can be represented as
$$
\vec{a}=  1\hat{x}+2\hat{y}+3\hat{z}
$$
where  $\hat{x}$, $\hat{y}$, and $\hat{z}$ are the orthonormal basis vectors (unit vectors along x-y-z directions).

However, these two things are different in one perspective.Among $x$,$y$,$z$, there is no concept as 'distance'. For example, in Fourier transform. we can say that the distance between $\omega_1 = 1$ and $\omega_2=2$ is smaller than the distance between $\omega_2$ and $\omega_3 =4$.

This algebra structure in Fourier transform definitely has some extra meaning, but what is it? Does it mean that the Fourier transform makes "more sense" compared with other orthonormal basis analyses where there is no easily-definable distance between basis functions?

Let $v(x, k)$ be an orthonormal set of basis functions such that
$$
\int_{-\infty}^{+\infty} v(x,a) v^*(x,b) dx = \delta_{a,b}
$$
where $\delta_{a,b}$ s the Kronecker delta.

I believe your question boils down to: for the Fourier Transform, we have
$$
v(x,a) = e^{jax}
$$
and we interpret $a$ as the frequency and $x$ as time; is there anything special about this orthonormal basis because we can interpret $a$ as a "distance" between these functions?

An example of another set of functions that are orthonormal (but over a finite interval) are the Legendre Polynomials:
$$
v(x,a) = P_a(x), \mbox{where} -1\le x \le +1
$$
where now $a$ is the order (degree) of the polynomial.

For these (and other polynomial bases), the "distance" is now the difference in polynomial degree.

Whether this means such functions are "less" special than the Fourier bases is a matter of opinion.


Generic response of an IIR filter from its transfer function How can I derive the general response of an IIR filter from its transfer function? I know that:

$$H(z)=\frac{1}{1+\sum\limits_{m=1}^N{a_m z^{-m}}}$$

Thus:

$$Y(z)=X(z)H(z)=\frac{X(z)}{1+\sum\limits_{m=1}^N{a_m z^{-m}}}$$

The general response is:

$$y[n]=-\sum\limits_{m=1}^N{a_m y[n-m]}+x[n]$$

Where the above expression comes from?

Thank you for your time.


Does using the Goertzel algorithm actually give better frequency resolution? I am reading this article, and I am getting a little confused by the author's liberal use of 'frequency resolution' regarding the Goertzel Algorithm.

Basic question: Does using the Goertzel algorithm actually give you more frequency resolution over a specific band of interest, or does it simply efficiently compute the FFT over only the specified band of interest, but at the same frequency resolution specified by sampling frequency divided by number of samples?

For example, lets say the $F_s$ is 100 KHz, (fixed) and the number of data samples $N$ is 10000. (Also fixed). If I compute a normal FFT, where FFT length is also $N$, my frequency resolution is $\frac{F_s}{N}$ as to be expected, and it will be equal to 10 Hz. This means my bins are separated by 10 Hz, from -50,000 Hz to 50,000 Hz.

Now let us say I want to use the Geortzel Algorithm to only look at frequencies in the range of say, 20,000-21,000 Hz. If I use the same $N$ for number of samples, and use the same $N$ for my FFT size, then what is my frequency resolution? Still 10 Hz? Or is it $\frac{21,000-20,000}{10000} = 0.1$Hz?

I have a feeling that I am not really increasing my frequency resolution, as much as simply interpolating points on the main lobe, by using the same $N$ in to evaluate the frequencies from 21,000 to 20,000 as I did from 0 to 50,000.

Is this a correct understanding?

Your understanding is correct.

Goertzel algorithms give almost exactly the same result as 1 bin of a DFT or FFT of the same length or number of samples (and where the FFT twiddle factors are generated by a trig recursion), when used for frequencies that are exactly integer periodic in the Goertzel length.  But many forms of the Goertzel algorithm provide only the magnitude, not the complex result or the phase of the FFT 1 bin result.  Also, computationally, a generic Goertzel may be numerically slightly less stable than a generic FFT.  For non-integer-periodic-in-aperture frequencies, the result is equivalent to a sinc interpolation between bins of a DFT or FFT of the same length (which can be slightly more accurate of an interpolation than the more typical parabolic interpolation of FFT results).

The interpolation might be said to increase resolution in a graphical sense (more plot points), or making it easier to visually spot the maxima, but not in the information theoretic sense, nor for better separating 2 close spaced spectral lines as 2 separate peaks.


What are real-valued and complex signals and why is the Fourier transform of a real-valued signal Hermitian? I have confusions in these concepts:

What are real-valued signals?
What is the difference between it and complex signals?
Why is it so that for real-valued signals, the spectrum of negative frequencies is the complex conjugate of the positive frequencies?

When one analyses data, outcomes both depends on the nature of the data, and their handling (eg a transformation, like the Fourier transform). One can consider binary data, categorical data, data in two dimensions, etc. To develop a theory and the related tools, it is useful to define where the d ata dwell.

A real signal is simply a signal that takes values that are real. The signal can be discrete or continuous, mono- or mulitdimensional, the observed values are real. The set of real values $\mathbb{R}$ is quite interesting. Since it is a continuous field, it allow a lot of "signal-like" mathematical operations: addition, product, convolution, and a lot of properties related to limits. $\mathbb{R}$ is one of the simplest mathematical structure allowing such operations.

But a natural tool to study signals, spectra,  is related to Fourier analysis. And Fourier uses complex arguments. Why shall we use complex transforms to study real signals? That is a complicated story, but let us say this is mandatory, to cope with linear and time-invariant systems. See for instance: Why cosine is not an eigen signal?

So one can also use complex-valued data. The complex field $\mathbb{C}$ is, maybe, the second simplest mathematical structure allowing the same operations as $\mathbb{R}$, without pain. And a real signal is just a subset of complex-valued data: take a complex object, whose imaginary part of zero, and you have a real signal.

For the last question: real signals are just a subset of complex signals. So when you analyze them (the reals) with a fully complex transformation, you can expect that the outcome will be restricted too. One of the restriction is that real signals have certain symmetries in the Fourier domain, that fully complex signals don't have


Why does an ideal high pass filter on this image generate gray values As seen from the image below, an ideal high pass filter was applied. Hence, stark changes in pixels from lighter to darker or vice versa show up on the filtered image. However, slower changes such as on the cheekbone appear to be this gray color.

What is this gray color? It seems to be quite prevalent in all my grayscale images.

Extreme positive values are black and extreme negative values are white. Values near zero are grey. Where there is not much change the high pass filter outputs values near zero. Hence the grey regions in the image.


What is the Use of the Opening Operator Just wondering what the opening operator is in image processing? What does it do? Just deletes bits of the image that you don't want?

Opening is the morphological operation used for example in removal of small particles or some noise in binary image. Additionally it is widely used in hand-writing recognition where you want only the 'skeleton' of your letter.

Mathematically it is a dilation of the erosion. Erosion sets current pixel to 0 if any pixel in neighbourhood is 0. Dilation on the other hand sets pixel to 1 if any neighbour is 1.

You can think of that in following way: first you do the erosion ($N$ steps), which is shrinking boundaries of your objects. After that you perform dilation (also $N$ steps), which is expanding boundaries of your objects. Because small ones were removed in erosion step you will bring back biggest elements.

Here is some image showing this operation:


Pattern recognition and Machine Learning reference My main reference for machine learning and pattern recognition is Pattern Recognition and Machine Learning by Bishop. But I want a complementary reference to this which has more exercises and explains the mathematical stuff in more details, e.g. formulas and proofs.

Another good book is:
 - Pattern Recognition by Theodoridis and Koutroumbas

It's about classifiers, features and clustering in the context of pattern recognition.

It contains all the mathematical stuff you want, and backs everything up with references to papers, books, etc. It has no exercises, so for that you need to look somewhere else.

Also see this StackOverflow question on the matter.


How to detect defects at a contour? I am detecting a contour of a smooth metal. I do binarization first then use openCV contour function. But sometimes the metal has defects such as dots or dusts at edge. This leads a unsmooth edge. Is there any way I can detect these defects with openCV function or some math calculation?
This is the contour result.

This is zoomed images for defected contour result.

Couple of approaches come to mind: you could pick an arbitrary "first" red pixel, then take a square of say 5×5 pixels around it, and simply figure out which quadratic function is a best fit for these, and infer the curvature from that.
Pick an arbitrary direction to start from that point and pick the next pixels; calculate the curvature of the 5×5 square surrounding that pixel. Find the absolute difference between last and current curvature; if it changed too much, you've found an abrupt corner.
In your specific example, applying to the whole image a 5×5 moving average filter would also suffice: Far as I can tell, there shouldn't be more than a relatively low number of red pixels in every 5×5 square, so if there's one with more, you've got your type of defect.
Other options include morphological classical tricks like erosion / dilation: use one operation to surround every red pixels with a fixed-diameter disk of red. Then do the same, with a diameter one linewidth (or more) larger disk, with the white in the picture. Check whether anything remains red afterwards.
PS: JPEG is the wrong image format / compression for line-style graphics. You'd want to work with PNG.


How to perform a block-wise frequency down-shift? I'm writing a Python function that shifts the frequency of an audio signal down to a specified range. I perform the downshift by multiplying the time series with a complex exponential, and then low-pass filter the result to remove the cyclic effect of the low frequencies shifting into the higher ones.

The function will be used to process an audio signal block by block in a loop. My question is, what will be the effect of starting the index $n$ of the complex exponetial at zero for each block?

$$ X \big(e^{j2 \pi (f-f_0)}\big) = \mathscr{F} \Big\{ x[n] \cdot e^{j2 \pi f_0 n} \Big\} $$

Where $n \in \{0,..N-1\}$, and $N$ is the block length.

Intuitively, I'm thinking there would be an inconsistent phase-shift between blocks, but I'm not sure.

is there a way to avoid this, without continually incrementing the index throughout the life of the program?

Here's the code:

from numpy import fft, pi, r_, exp, int16, array, fromstring, zeros
from scipy.signal import firwin, lfilter

def shiftFreq(data, cutoff, fs, hist=None):
    """ Downshift audio data

    Parameters
    ------------
    data: str
        Binary data string of 16-bit audio data
    cutoff: float
        The new max frequency component after downshift and filtering
    fs: float
        Sample rate of audio data
    hist: ndarray
        Filter delay line

    """
    data = fromstring(data, dtype=int16)
    N = len(data)
    num_taps = 128
    delay = num_taps - 2
    nyquist = fs/2

    # initialize history
    if hist is None:
        hist = zeros(delay, dtype=data.dtype)

    # Down-shift the data
    shift = nyquist - cutoff
    n = r_[0: N]
    new_data = (data * exp(-1j * 2 * pi * shift * n/fs)).real

    # Apply lowpass filter to remove cyclic effect
    fir = firwin(num_taps, cutoff, nyq=nyquist)
    new_data, hist = lfilter(fir, 1, new_data, zi=hist)

    return int16(new_data).tostring(), hist

As an aside, are there other caveats I should be aware when performing a down-shift, that I haven't accounted for here?

As an aside, are there other caveats I should be aware when performing
  a down-shift, that I haven't accounted for here?

Yes there's one and let me put one of the most important caveats when performing frequency modifications. An acoustic sound signal (in the form of a music, speech etc) will physically exist in a so called harmonic structure (for more indepth discussions see sound synthesis, sound analysis etc.) whose immediate implication is that the perception of the frequency content is in octaves rather than linear frequency decades (an octave is the span of frequencies between $f_0$ and $2 f_0$).

Now, eventhough for most technical applications (such as in those communication algorithms) a band of frequencies can be shifted to any desired frequency without a loss of information (which forms the essence of the modulation theory) when it comes to audio perception this does not hold.

It's easy to see that the harmonic structure of the audio band will be distorted when a complete band is shifted to other frequencies in a linear manner. As a solid example consider a primitive audio waveform composed of $440$ Hz, $880$ Hz and $1320$ Hz (these are $f_0$, $2f_0$ and $3f_0$ of an absurd instrument, more natural instruments will have more complicated structures including inharmonicities etc.)

Now if you linearly shift this waveform base from $f_0$ to $1.5 f_0$ then the new frequency distribution will be $1.5f_0$ , $2.5f_0$ and $3.5f_0$ which are no more harmonics. This creates an important problem in natural pitch shifting applications.


Frequency response of Low pass filter and High pass filter The transfer function of a Low pass filter is H(w). From this I want to develop a high Pass filter. I read in DSP by Proakis, Sec 4.5, that the high pass filter can be obtained by translating H(w) by  PI radians. I did not get how it will form a High Pass Filter? I think tt should form a Bandpass

In previous sections of the book, the fact that a discrete-time signal's spectrum is periodic may have been mentioned. It can be described formally as follows:
$$X(e^{j\omega})=\frac1{T} \sum_{k=-\infty}^{\infty}X_C\biggr(j\biggr(\frac{\omega}{T}-\frac{2\pi k}{T}\biggr)\biggr)$$
being $X_C(j\omega)$ the Fourier Transform of the continuous signal, and $T=1/f_s$ the sampling period. That expression is actually saying that the discrete-time spectrum can be computed from the continuous spectrum just by expanding its width $T$ times, and repeating this every $\frac{2\pi}{T}$. To keep things simple, we can normalize the frequency so that every signal's spectrum is periodic every $2\pi$, no matter how we sampled it.

Keeping this in mind, a low pass filter isn't just a "step" around $\omega=0$, but also around $\omega=2\pi,4\pi,...$ and, of course, their negative counterparts. I mean low frequencies repeat around $\omega=0, 2\pi,4\pi,...$, while high frequencies do around $\omega=-\pi,\pi,3\pi,...$

Shifting its frequency response by $\pi$, we get $H(e^{j(\omega-\pi)})$, which would look like the following:

Remember that $\omega=\pi$ corresponds to your Nyquist frequency (= half the sampling frequency).

Image credit: Oppenheim & Schafer, Discrete-time signal processing


Is the zero part of a RRC filter generated signal sent? When you use a RRC filter you add extra symbols at the end of your generated data because of the delay of the filter. When the data passes through the filter it gets up-sampled and the delay is introduced. However, the delay also introduces some zero, or near zero, valued parts to the signal, see figure. My doubt was, at the moment of sending the signal, are this parts also sent? Or you only take the non zero value parts?

That depends. But typically, yes, you'll flush out your transmit signal; that basically "costs" nearly nothing, but makes sure your spectral mask stays clean.

At the receiver, this part still contributes to the amplitude of the matched filter output for the last couple of symbols, so, albeit not contributing very much, omitting it would just introduce an error to these symbols.


Auto correlation correspondence to original time domain window I'm trying to figure out how exactly an auto correlation corresponds to a time domain signal.

Now I'm trying to find pitch periods in an audio file and 99% of the time I'm getting it spot on.  Unfortunately however that 1% of the time is throwing me some annoying issues.

Now at present I scan through my window for the highest peak (ie the highest abs( sample value )).  This seems to match up really nicely in most cases.

Now when I auto correlate I would expect the central peak to be this peak that I have identified (It will always have a positive number if i read it right).  Is this correct?

Secondly, A peak that is n samples away from the central peak of the auto correlation will be n samples away from the peak I've identified in time domain, ok?

Sooo, if this is correct then I'm having a strange problem.  I'm finding a peak in the auto correlation 91 samples from the center of my auto correlation.  Unfortunately when I look at the signal in the time domain (inside an audio editing app) then my central peak corresponds perfectly to a negative peak in the audio.  However, there are no peaks 91 samples away from that peak.

If I look forward the next real peak is 142 samples away and the previous peak to it is 75 samples away.  I do notice a very tiny peak 91 samples or so in front of my peak but the sample never even goes negative.  So why does this end up showing up as by far the strongest peak in the auto correlation?

Any help with this would be hugely appreciated!

Autocorrelation is not about finding the distance between individual peaks.  It is more about finding those lag distances that minimize the averaged squared delta between everything, all the peaks, all the valleys, all the flat spots, all in combination, and etc.  Because of this averaging over the entire window, the lag distance may not correspond to the distance between any single pair of point phenomena.


Calculating an HSV histogram I have to compute an HSV histogram of a colored image.  Is it same as an H-S histogram?

I have found H-S histogram computation code everywhere. What is the method to change an H-S histogram to an HSV?

Since you want an histogram, I believe you want to plot it. In GNU Octave you may use hist3 and make three histograms (one for each pair: HxS, HxV and SxV). The code follows bellow:

pkg load statistics
I = imread ('lena.png');
sI = size (I);
Ihsv = rgb2hsv (I);
figure ("position",get(0,"screensize")./[1 1 2 3]);
subplot (1,3,1);
hist3 ( reshape (permute (Ihsv(:,:,[1,2]), [3, 1, 2]), 2, prod(sI(1:2)))' );
title ('H vs S'); xlabel ('H'); ylabel ('S');
subplot (1,3,2);
hist3 ( reshape (permute (Ihsv(:,:,[1,3]), [3, 1, 2]), 2, prod(sI(1:2)))' );
title ('H vs V'); xlabel ('H'); ylabel ('V');
subplot (1,3,3);
hist3 ( reshape (permute (Ihsv(:,:,[2,3]), [3, 1, 2]), 2, prod(sI(1:2)))' );
title ('S vs V'); xlabel ('S'); ylabel ('V');


Plotting antenna pattern in MATLAB I have a beamformed output, say

y = weights * signal,

where the signal is gaussian,i.e signal = randn( 1 , 1e5 ).

Now I want to visualise the output in angular domain, i.e I want to plot the gain of the resulting signal y as a function of elevation angles. I am considering the mainlobe of 'signal' to be pointing in (0 azimuth,0elevation). The weights are combination of amplitude weights and phase weights(exp(phase terms)). I am computing the antenna pattern in the following way:

pattern = abs(arrayFactor.* y).^2;

Is it correct to plot the antenna pattern of such a signal y by multiplying it by the array factor and then plotting its squared magnitude?(in case of isotopic element pattern.)

I would highly appreciate any input.

Denote the array manifold (what you call array pattern) by a $N \times 1$ vector ${\bf{a}}\left ( \theta \right )$ where $\theta$ is the DOA and $N$ is the number of array elements. Further denote by the $N \times 1$ vector $\bf{x}_k$ the $k$-th snapshot of the array, then you are measuring the following signal:

\begin{equation}
{\bf{x}_k}= {\bf{a}}\left ( \theta \right ) s_k + {\bf{n}}_k
\end{equation}

where $s_k$ is the signal reflecting from the target and $n$ is a random noise.

Denote the beamformer by a $N \times 1 $ vector of weights ${\bf{w}}$, then the array response after beamforming would be

\begin{equation}
r(\theta) = {\bf{w}}^H {\bf{x_k}},
\end{equation}

Note that $r$ is a scalar for a given DOA $\theta$ and $\bf{w}$ is a function of $\theta$. Now, to get the array response for a variety of DOAs just change $\theta$.


Energy per bit ($E_b$) in BPSK with non-unity channel gain I would like to know what is exact and correct definition of $E_b$ (energy per bit) in digital communication. Consider I have a model such as      r=a*x+n where a is the channel gain, x is a (+1 or -1) BPSK symbol and n is zero mean AWGN noise. I am defining $E_b$ to be (assuming a=10 and x given:

a=10;
Eb= sum((a*x).^2)/length(x);

However this method seems to be failing since $\sigma$ can be related to $E_b$ using $\sigma=\sqrt(E_b/(2*\text{SNR}))$. So if a is increased the noise power is increased as well and this should not be the case. Therefore should the definition of $E_b$ be like following:

   a=10;
   Eb= sum((x).^2)/length(x);

Meaning that is $E_b$ normalized to the a (channel gain). In this case then $E_b$ is 1 and as a increases the noise power won't be increasing. Can anyone elaborate on this? Should the noise power increases as the amplitude of signal increases? Then in this case performance gets worse and worse as signal amplitude is increased. Is this a correct statement?

The energy per bit $E_b$ means one thing and one thing only: how many joules is the transmitter spending, on average, per information bit transmitted. Note that the channel has absolutely nothing to do with $E_b$.

However, your question involves the SNR too. In a channel with non-unit gain, you need to decide where you'll measure the SNR.

You can define the SNR at the transmitter: $\text{SNR}=2E_b/N_0$ for BPSK. Obviously, the BER depends on both the SNR and the channel gain.

You can also define the SNR at the matched filter's output: $\text{SNR}=2G^2E_b/N_0$ for BPSK, for a channel with gain $G$. In this case, the BER depends exclusively on the SNR.

When $G$ is random with known probability density, it is common to average the SNR over $G$, in order to obtain the expected BER over many channel realizations.


Subsampling data: terminology I'm writing software and the data reporting frequency is slower than the acquisition cycle, so effectively sub-sampling is being done.  Different data types require different methods of sub-sampling which I'm tempted to call a Gating method where oversampled data in one time domain goes through a Gate to another time domain. Windowing also comes to mind.

I hesitate though as although these terms imply an interval of focus/ interest, they might not necessarily imply continuity between intervals; back to back intervals where no input is ignored.

Such methods might be:

Average/Median since last report (or low pass filter)
Min/Max since last report
Sticky flags (a bit field who's bits are cleared when read)

Is there a better term than Gate or Window?

Note 1: This is borderline a programming question or code review, but I think still more for signal processing.

Note 2: Although sub-sampling causes signal artifacts such as aliasing.  This is not part of the question.

This is more likely a question of what the data reporting application is trying to achieve, rather than a DSP question.

Yes, there is a different term for Gate or Window. In fact, if you "subsample" and low pass filter, as per your first intuition, then you are effectively doing decimation. For more information please see: http://en.wikipedia.org/wiki/Decimation_%28signal_processing%29

The rest of your suggestions would be equivalent to subsampling followed by non-linear filtering such as median, min, max which would have the corresponding effect on the reported data.

For example, a sliding window median filter would "absorb" any spikes in the reported data and the min and max filters would tend to report the highest or lowest values. For more information please see: http://en.wikipedia.org/wiki/Nonlinear_filter and http://en.wikipedia.org/wiki/Median_filter

So, it is really a matter of what does that GUI element try to convey to its user.

Hope this helps.

Supplementary Information:
Yes decimation does imply a stable sampling frequency but most DSP operations can be defined over irregular domains as well.

Irrespectively of that however, I think that we are indirectly looking at the real problem here which is one of visualisation.

The real problem, or at least my understanding of it so far, is that you are trying to convey a comprehensive view of a window of the history of a signal at a rate lower than its actual rate of change. Furthermore, this signal is not used in any further computation. If that is mostly accurate and assuming that your signal is coming from a stochastic process then no single value will be able to represent this history reliably.

Therefore, I would suggest two alternatives:

To use a single boxplot that would represent your last signal's frame
history in a comprehensive way including its min, max, median and
quartile information. For more information please see:
http://en.wikipedia.org/wiki/Box_plot
To convey the signal's history to your user via a Sparkline which
could furthermore be annotated with min,max,median lines. Please note that Sparklines can also be created via the use of special fonts, which makes them really easy to quickly insert to a GUI. For more
information please see: http://en.wikipedia.org/wiki/Sparkline and https://www.fontfont.com/how-to-use-ff-chartwell or http://nsfmc.github.io/chartjunk/

Hope this helps


MATLAB fvtool: Setting the passband ripple below 0dB for a FIR Filter I'm new to using filter tools in MATLAB and I'm having a bit of difficulty. I'm designing a lowpass FIR filter and I'm trying to keep the passband ripple below 0dB in the magnitude response. So far, the passband ripple travels up above 0dB the same amount that it travels below. Any ideas on how to accomplish this using fvtool? I've been reading up on fvtool all day and I still can't figure it out. Thanks in advance for your help. Also, I apologize for this being so wordy for such a simple questions.

From your question I conclude that it is an equi-ripple filter that you designed, probably using the Parks-McClellan algorithm. Such an equi-ripple FIR filter has a magnitude response alternating between the values $1+\delta_p$ and $1-\delta_p$ in the pass band (where $\delta_p$ is the maximum approximation error in the pass band), and in the stop band it also has ripples of equal height $\delta_s$, where $\delta_s$ is the maximum stop band error. Note that all quantities I use here are linear, not in dB.

So, as mentioned by SleuthEye in a comment, after designing such a filter you can simply scale it to get the desired behavior. If $h[n]$ are your original filter coefficients, your desired filter has the coefficients

$$\tilde{h}[n]=\frac{1}{1+\delta_p}h[n]\tag{1}$$

The filter in $(1)$ will have a pass band magnitude alternating between the values $1$ (i.e., 0 dB) and $(1-\delta_p)/(1+\delta_p)$, and its minimum stopband attenuation will be $\delta_s/(1+\delta_p)$.


Formula to compute output of a linear system The impulse response is $h(t)=e^{-t}u(t)$  and input signal is $x(t)=1 + \frac{1}{2}cos(400πt)$

I want to find y(t)..

I tried to use $\int^{\infty}_{-\infty} h(τ)x(t-τ)dτ$ ,however i know that   $\int^{\infty}_{-\infty} h(t-τ)x(τ)dτ$  is equal.
Is this the correct formula?
I tried to solve it but couldnt find come up to a solution..

u(t) is the http://en.wikipedia.org/wiki/Heaviside_step_function

It's correct. Here's the general way how you can solve such a problem – in this specific example it would actually be easier to just solve the integral, but that seems indeed like homework you should be doing yourself.

I'll call the step function $\Theta$, that's how we always call it in physics – hope it doesn't confuse, and the convolution
$$
  (f\star g)(a) := \int_{\{(t,t')\colon t+t'=a\}}\!\!\!\!\!\!\!\mathrm{d}(t,t')\;f(t)\cdot g(t') = \int_\mathbb{R}\!\mathrm{d}\tau\; f(\tau)\cdot g(t-\tau)
$$
which is just the definition that you use. We know the convolution theorem
$$
  \mathrm{FT}(f\star g)(\omega) = \frac1{\sqrt{2\pi}}\mathrm{FT}(f)(\omega)\cdot\mathrm{FT}(g)(\omega)
$$
(the factor $\frac1{\sqrt{2\pi}}$ depends on your definition of the Fourier transform). To apply this to your problem first calculate
$$
  \mathrm{FT}(h)(\omega) = \frac1{\sqrt{2\pi}}\int_\mathbb{R}\!\mathrm{d}t\; e^{-t}\Theta(t)\cdot e^{-i\omega t}
 = \frac1{\sqrt{2\pi}}\int\limits_0^\infty\!\mathrm{d}t\; e^{-t}\cdot e^{-i\omega t}
$$
since the Heaviside function cancels everthing for $t<0$, then this is
$$
 = \frac1{\sqrt{2\pi}}\int\limits_0^\infty\!\mathrm{d}t\; e^{-(1+i\omega) t}
 = \frac1{\sqrt{2\pi}}\left.\left[\frac{e^{-(1+i\omega) t}}{-(1+i\omega)}\right]\!\!\right|_{t=0}^{\infty}
 = \frac1{\sqrt{2\pi}}\Bigl(\frac{0}1 - \frac{e^0}{-(1+i\omega)}\Bigr)
$$
which is just $=1/(\sqrt{2\pi}\cdot(1+i\omega))$. Meanwhile,
$$
  \mathrm{FT}(x)(\omega) = \mathrm{FT}(\lambda t.1)(\omega) + \tfrac12\mathrm{FT}(\lambda t.\cos(400\pi t))(\omega).
$$
These Fourier transforms only exist in a distributional sense; it is known (can easily be shown by inverse $\mathrm{FT}$) that
$$
  \mathrm{FT}(\lambda t.1) = \sqrt{2\pi}\;\delta
$$
and
$$
  \mathrm{FT}(\cos)(\omega) = \frac{\sqrt{2\pi}}{2}\bigl(\delta(\omega-1)+\delta(\omega+1)\bigr).
$$
So we have
$$\begin{aligned}
  \mathrm{FT}(y)(t) =& \frac{\sqrt{2\pi}\bigl(\delta(\omega)+\tfrac14\delta(\omega/400\pi-1)+\tfrac14\delta(\omega/400\pi+1)\bigr)}{\sqrt{2\pi}(1+i\omega)}
  \\=& \frac{\bigl(\delta(\omega)+\tfrac14\delta(\omega/400\pi-1)+\tfrac14\delta(\omega/400\pi+1)\bigr)}{1+i\omega}
\end{aligned}$$
and therefore
$$\begin{aligned}
  y(t) =& \mathrm{IFT}(\ldots)(t)
   = \frac1{\sqrt{2\pi}}\int_\mathbb{R}\!\mathrm{d}\omega\; e^{i\omega t}\frac{\bigl(\delta(\omega)+\tfrac14\delta(\omega/400\pi-1)+\tfrac14\delta(\omega/400\pi+1)\bigr)}{1+i\omega}
 \\=& \frac{e^0}{\sqrt{2\pi}\cdot 1} + \frac{e^{-400 i\pi t}}{\sqrt{2\pi}(1+\tfrac1{400\pi}i)} + \frac{e^{400 i\pi t}}{\sqrt{2\pi}(1-\tfrac1{400\pi}i)}
\end{aligned}$$
We could know calculate this as an exact expression of trigonometric functions by writing out the $1\pm\tfrac{i}{400\pi}$ as a polar complex number.


Probability of random numbers This is a question I found on internet under DSP section so that's why I am posting it here. Help me understand it please.

  A computer adds 1000 random numbers that have each been rounded off to the nearest 10$^{\rm th}$. Find the probability that the total round-off error for the sum is $\ge 1$.

How probability of any system can be greater than 1?

The total roundoff error for the sum of $N$ numbers is:

$$
S = \sum_{i=0}^{N-1} E_i
$$

The roundoff error for the $i$-th number is represented by the random variable $E_i$. If we assume that the random number generator used by the computer yields numbers $X_i$ taken from a uniform distribution, then the difference between each $X_i$ and the nearest tenth (which is the roundoff error $E_i$) is uniformly distributed on the interval $(-\frac{0.1}{2}, \frac{0.1}{2}) = (-0.05, 0.05)$.

What we're concerned with, though, is the distribution of $S$. Since $S$ is the sum of $N$ independent, identically distributed (iid) random variables, then via the central limit theorem, as $N \to \infty$, $S$ will tend to a Gaussian distribution. If we assume that your case of $N=1000$ is "large enough" for the Gaussian assumption to hold, we can easily estimate the probability that you seek. It's certainly possible to exactly calculate the distribution of $S$, but the Gaussian assumption is likely close enough for most applications with such large $N$.

A Gaussian distribution is characterized by its first two moments, so if we can find those for $S$, then we have all the information we need. These are easy to calculate for a sum of iid random variables. The mean of $S$ is equal to:

$$
\mathbb{E}(S) = \sum_{i=0}^{N-1} \mathbb{E}(E_i) = 0
$$

The variance of $S$ is equal to:

$$
\mathbb{E}\left((S - \mathbb{E}(S))^2\right) = \sum_{i=0}^{N-1} \mathbb{E}\left((E_i - \mathbb{E}(E_i))^2\right)
$$

Recall that the random variables $E_i$ are distributed uniformly. It is well known that the uniform distribution over the interval $(a,b)$ has variance $\frac{1}{12}(b-a)^2$. For this case, that yields a variance $\sigma_{E_i}^2 = \frac{0.01}{12}$. Therefore, the variance of the total roundoff error $S$ is $\sigma_{S}^2 = \frac{0.01N}{12}$.

So in summary, we can approximate $S$'s distribution as Gaussian with mean zero and variance $\sigma_{S}^2 = \frac{0.01N}{12}$. Based on those parameters, you can easily calculate the estimated probability distribution function (pdf), then integrate that result to arrive at whatever probability you seek. The probability that there is a total roundoff error with magnitude greater than one would be:

$$
\begin{align}
P(|S| > 1) &= P(S>1 \lor S < -1) \\
&= 1 - P(-1 < S < 1) \\
&= 1 - \int_{-1}^{1}f_S(s)ds
\end{align}
$$

where $f_S(s)$ is the Gaussian distribution's pdf that we arrived at before.


Image Geometry Terms I getting confused over some basics of image geometry terms due to different definitions, implementations and different outputs in the internet, here is an example. The terms are:

Image derivative vs Image gradient.

Image directional derivative/gradient.

Image Laplacian  vs   image second derivative   vs image divergence.

Image curvature (how it is related to the other terms?).

Image normal vs image gradient magnitude.

it would be nice if someone explains these terms (or provide a reliable source so I can refer for more questions) in a simple way providing:

Simple math symbols.
Meaning of these terms (what can tell us about the input image?).
Simple implementation e.g. simple matlab (without buildin functions) or c++ code.
An input and output images.

All the best!

Ibraheem

Image derivative vs Image gradient.

Same thing.

It tells you how much local contrast there is between adjacent regions in an image. The image gradient gives you the slope of the contrast towards its greatest magnitude.

  Image directional derivative/gradient.

Same thing. In a simple $f(x)$, direction is implied. The derivative is always taken with respect to $x$ and the $x$ axis points to the right (if we are talking about a spatial dimension). In an $f(x,y)$ or more generally, any $f$ of more than one variables, the derivative can be obtained towards (or parallel to) a direction. If you don't specify the direction of differentiation, you get the  derivative (as above) which tells you the direction of higher slope.

  Image Laplacian vs image second derivative vs image divergence.

You need to understand Divergence. Divergence tells you how much "outgoing flow" there is in a vector field. Imagine a fountain, the water comes up through a pipe in a small pool. "New" water comes in to the pool and pushes the existing water out of its way. Assuming that the pool is perfectly levelled, "new" water is pushing existing water outwards without preference to a specific direction.Isotropically. If we were to represent the velocity field as a vector field, at the point of the fountain we would see vectors going out, radially, towards all directions, like a little star. The Divergence of that point (the source) would have positive maximum value. If the exact opposite happened (i.e. the water was flowing IN a sink), then the divergence at that point would have negative maximum value.

So, how do you estimate the Divergence at some point? You take a "circle" (or sphere or hypersphere in higher dimensions), you position it at a particular point in space and then you estimate how much flow goes in and out through its boundary. Why circle / sphere / hypersphere? Because Isotropically. How do you do this "boundary" thing? With a surface integral. Where do you get the "flow" from? From the vector field.

The Laplacian is the Divergence of the gradient. Remember that Divergence expects a vector field at its input.

The Laplacian of an image depends to an extent on what the image depicts. If we are talking about simple visible light images, the Laplacian gives you areas that either reflect light or absorb light. It can be used to find specular reflections for example or spherical regions that bounce light towards all directions (concave or convex).

The second derivative is the second derivative (please see earlier answers). It tells you the slope of the slope of "brightnes" and it can be used to find local maxima and local minima.

  Image curvature (how it is related to the other terms?).

Which image curvature? You might have curvature of features ON the image (e.g. the curvature of a fence).

Given the discussion so far, you probably mean "curvature" as it emerges in the third dimension (reflectivity or brightness). In that case, curvature is more related to the Laplacian (but it could also be related to gradient towards a direction, anything that implies curvature).

  Image normal vs image gradient magnitude.

A Normal usually refers to a vector. A vector has direction and magnitude.

The image normal is a vector that is perpendicular to some notion of curvature of the image.

Usually, normals are made to have a magnitude of one because we are only interested in the direction they are pointing towards. But that is not necessarily the case with any normal.

The image gradient tells you the highest slope departing or arriving at a given point in the image. Therefore, the image gradient is a vector field which could be used to derive the normal vector field. But if you simply take the magnitude of it, all that you have is an estimate of local contrast towards the direction of maximum difference.

My advice would be for you to revisit vector calculus and linear algebra because these are the fundamental concepts that image processing is re-using. Divergence is not redefined in image processing. This will also help you in terms of "code". All major scientific computing platforms have these operators implemented, it is just a matter of re-using them rather than re-writing them.

The other thing to note is that imaging has context. Divergence over visible light photography is different in interpertation than infrared photography or X-Ray photography or other types of imaging. In visible photography there are certain assumptions that are often made to simplify the mathematics (e.g. in the way that light reflects off of objects). These assumptions make it "easy" to use something like Divergence to estimate curvature but there are still limitations. (And you might end up requiring more than one images to truly estimate "curvature"). In visible light photography, bright / dark, usually means light/shadow, but in infrared photography, bright/dark could also mean near/far.

Hope this helps.