Signal Processing for Discrete Events This is less a specific question, and more... "what is this called" and "where can I read more about it"?

Several times in my career, I've had to work with datasets with "discrete events". Lets say customers passing an infrared sensor at a store front.

Normally, signal processing would be concerned with a continuous time signal, like the voltage on the IR receiver... but my question pertains to the event of someone walking through the door. The "signal" (in this case) is more a list of times when someone walked through the gate.

For example:

# Unix Timestamps
entry_times = [1490648032, 1490648102, 1490648591, ..., ...]

The "crude" way to treat this signal is to bin the entry times to hours of the day and interpolate that dataset for a rough idea of customers vs time. However, I think that there's more to this problem.

Entry and exit are the same type of signal, but could DSP help differentiate? (I think not, but I don't know!)
Perhaps we want to count a mother and her two children as a single "customer". Assuming we get three events in a five second window, can we somehow group that together?

If we count the number of entries and exit per hour, this gets us a rough graph of customer throughput, and we can interpolate that to get more resolution... but how could we create that on the fly?

What is the name of this type of signal? Are there any signal processing algorithms or techniques that would help answer my questions?

To do something about the "Mother with kids"-Problem, you could implement an algorithm that sets a value to 1 or 0 based on whether the sensor was activated in a timeframe of ~5 seconds. Saving this value every 5 seconds leads to equidistant samples in time domain.

The coming in/going out-problem could be solved using two IR sensors that customers have to pass one after another to enter or leave. By checking whether the activation time difference is positive or negative, you can find out if they were coming or leaving.

I would not recommend interpolation as customer appearance is not necessary bound to any pattern or formula but somewhat random since it depends on a lot of things not easy to predict.
In this case, consider using histograms for fixed time intervals (15/30/60 minutes) to visualize customer spread (rush hours ...) .

Assuming you take aforementioned samples every 5 seconds, you are working with value- and time discrete signals which may or may not be deterministic/periodic (predictable) but will be random/stochastic to some extend.


What is the spectrum of modulated noise? I use a white-noise generator in an rf bridge to measure antenna impedance.
The receiver is an ordinary AM/CW/SSB communications type, nothing fancy, nothing 'coherent'.

As I approach the null of the bridge the signal is masked by receiver noise.  To make it stand out I 'modulate' the noise with a 1kHz square wave (the noise source is switched on/off at about 1kHz.)  I 'null' on the 1kHz audio, not the received noise power.

The question:  What is the spectrum of the modulated noise and how is it different from that of the 'raw' noise?

Short answer: most likely white noise.

Long answer:

Amplitude modulation using a sinusoidal wave shifts the original spectrum in frequency by the sinusoids frequency.  For instance, if there is a source with a square spectrum spanning -100Hz to 100Hz, and it is modulated by a 1000Hz sinusoid, it will have shifted the spectrum in both positive and negative frequency.  Thew new spectrums will show two sources of energy from -1100Hz to -900Hz and 900Hz to 1100Hz.

To answer your question, we also should consider that their source signal is white-noise, and the carrier signal is a square wave rather than a sinusoid.  A square wave, rather than having a single peak in the frequency domain (such as a sinusoid does), has multiple harmonics of descending strength.  Modulating with a square wave will create multiple images of your source signal to appear, each centered around the square wave's harmonics and scaled by the half the strength of the harmonic.

If your input signal is true white noise (really only observed mathematically), the shifted noise signals would only add up to more white noise.  But, if you have band-limited noise, then you will see copies of the noise spectrum shifted to the square wave's harmonics.


How do I detect the sharp edges in the following image? I am looking to count the sharp pointy triangular things in the image. However I was unable to do so. One approach I did try was take a convex hull of the contour and then find convex defects but it led to lot of irregularities. You can find the it in the labelled image below.

Ok I solved this problem using a simple solution. After I computed my convex defects in my countours, I performed a point polygon test and found the distance between the point in the convex defect and the hull. And I gave a threshold after looking at the distance of the points and all the sharp points were in the threshold greater than approximately 30. So this filtered my points


Path to becoming an audio engineer Can someone write down steps an totally beginner in EE and DSP needs to pass to become a audio engineer? From which subject to start with (math, electronics), and with which subjects to continue this road? My goal is to get a good knowledge in audio signal processing, things like creating sounds (sound synthesis), sound transformation, mixdown audio, understand fundamental things. I will be very grateful if someone give me instruction from where to start and what path to follow, and on what to mostly pay attention, I have a great desire to learn all of that, I have 6 years experience in electronic music production. Thank you very much!

(TLDR: if you want to dive in I might recommend Will Pirkle's books Designing Audio Effect Plug-Ins in C++ and Designing Software Synthesizer Plug-Ins in C++ but you might be completely overwhelmed depending on your background)

First off, if you want to create new sounds, I love Max/MSP and that will get you creating new sounds quickly but it will not teach you DSP.

Before a mod closes this as off-topic or too open-ended, I will say this:

There is no recipe, no set list of steps. Period.

Or as Umberto Eco wrote in his famous thriller "Foucault's Pendulum"....there is no map.

Take that for what you will. I'm sure everybody's story who works with audio is different. I grew up messing around with 4-track tape recorders, and since I grew up with computers, started playing with digital audio workstations like Audacity and ProTools (Ardour during my hardcore linux purist phase), making demos and playing in bands. During this time I also tried my hand at BASIC and later scripting languages. I got interested in theoretical physics in high school, and majored in physics in college (where I made a lot of demo recordings with bands and learned how to solve problems, but not a lot of domain specific knowledge (which I think is a good thing)).

I didn't want to be a particle physicist and had met a professor who worked in architectural acoustics. I found a grad program that was able to pay me to think about sound and space (mostly psychoacoustics) for 4 years and earn my MS and PhD.

At the end of my PhD I met (through a random LinkedIn message I initiated) the CTO of a startup who was willing to transfer my skills in psychoacoustics to programming microcontrollers to do effects and equalization. My job got me started down the path of really doing synthesis and processing. I had (almost) no idea what I was doing, and that's when I really learned. I'm just starting on this journey of being able to digitally manipulate sound, and I find that there's just so much out there and it is totally fascinating.

Every day I realize I know so little, but take in what I can each day (and try not to get frustrated that I don't move faster). Understanding is a slow process that cannot be expedited.

That being said, if I had to do it all over again knowing what I know now (and how could I have really?), I would have taken as much math and programming as I could have in high school (which wasn't offered in Catholic school) and majored in electrical engineering (which wasn't offered at the liberal arts college I attended). EE is really where all this stuff stems from. I'm biased but I would say it doesn't hurt to study something about hearing science and psychoacoustics too. Music appreciation can fuel your enthusiasm.


Frequency filtering with a DFT and meaning of removing complex conjugates I've seen many implementations to filter out frequency components of some time domain signal by performing a DFT, zeroing the unwanted frequency bins, and performing the IDFT to get the filtered signal, but more than a few times I've seen the results of the DFT treated as frequency components from 0 to the sampling rate instead of 0 to half the sampling rate then the complex conjugates. Take for example this bandpass implementation by a group at MIT: https://github.com/diego898/matlab-utils/blob/master/toolbox/EVM_Matlab/ideal_bandpassing.m

It's defining the frequency bins and corresponding bandpass filter as

Freq = 1:n;
Freq = (Freq-1)/n*samplingRate;
mask = Freq > wl & Freq < wh;

n being the number of samples and wl and wh being the low and high frequencies respectively.

My first intuition is that that's incorrect, but it makes me wonder if I'm missing something since I've seen it done this way multiple times.

But the more interesting question to me is, assuming it is wrong: what physical meaning, if any, does removing one but not both of these complex conjugates have? For example, if I used the above bandpass implementation to make a low-pass filter which filtered out one half the sampling frequency, I would keep the DC component and all the first complex values. But then I would remove all of their conjugates. While I don't actually think that's the low pass filter intended, it does still have the effect of smoothing out the signal, but how would you describe that? I don't imagine it's still considered a temporal filter in that case.

Here is a different way to look at it:

The violation of complex conjugate symmetry produces and "error". That error is, that the impulse response becomes complex. By applying the "real()" operator on the inverse FFT, the error is eliminated again.


Implementation of Wikipedia Equation for the DFT I was writing a simple fourier transform implementation and looked at the DFT equation on wikipedia for reference, when I noticed that I was doing something differently, and after thinking about it felt that the wikipedia version must be wrong because it's very simple to think of a signal that when fourier transformed (with that equation) will return an incorrect spectrum: Because the equation wraps the signal around the complex plane only once (due to the $n/N$ with $0<n<N-1$), any signal that is periodic an even number of times (while wrapping the complex plane) will have no spectrum as the usual peaks (while going around the unit circle) that would appear during a DFT will cancel each other out (when an even number of them appear).

To check this I wrote some code which produced the following image, which seems to confirm what my thoughts.

"Time using equation" uses the equation
$$ X_f = \sum_{n=0}^{N-1} x_n (\cos(2\pi f t_n) - i \sin(2\pi f t_n)) $$
with $t$ a vector of time (so the time $t_n$ at which $x_n$ was sampled for example). It can be found in the function ft below.

The wikipedia equation, linked above, is copied here for reference:
$$ X_f = \sum_{n=0}^{N-1} x_n \left(\cos\left(2\pi f \frac{n}{N}\right) - i\sin\left(2\pi f \frac{n}{N}\right)\right) $$
It can be found in the function ft2.

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t))
        compl = np.sum(- s * np.sin(2*np.pi*freq * t))
        tmpsum = (real**2 + compl**2) ** 0.5
        S.append(tmpsum)
    return S, freqs

def ft2(s, fs):  # Using wikipedia equation
    nump=len(s)
    freq_step = fs / nump
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for i, freq in enumerate(freqs):
        real = np.sum(s * np.cos(2*np.pi*freq * i/nump))
        compl = np.sum(- s * np.sin(2*np.pi*freq * i/nump))
        tmpsum = (real**2 + compl**2) ** 0.5
        S.append(tmpsum)
    return S, freqs

def main():
    f = 5
    fs = 100
    t = np.linspace(0, 2, 200)
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs)

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S, freqs = ft2(y, fs)
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()

Obviously it seems rather unlikely that I would have randomly found an error on such a high profile wiki page. But I can't see a mistake in what I've done?

You have a bug in ft2.  You are incrementing i, and freq together.  That's not how you want your summation to work.  I messed around with fixing it, but it got messy.  I decided to rewrite it from a discrete perspective instead of trying to use the continuous terminology.  In the DFT, the sampling rate is irrelevant.  What matters is how many samples are used (N).  The bin numbers (k) then correspond to frequency in units of cycles per frame.  I tried to keep you code as intact as possible so it would remain easily comprehensible to you.  I also unfurled the DFT calculation loops to hopefully reveal their nature a little bit better.

Hope this helps.

Ced

import numpy as np
import matplotlib.pyplot as plt

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t))
        compl = np.sum(- s * np.sin(2*np.pi*freq * t))
        tmpsum = (real**2 + compl**2) ** 0.5
        S.append(tmpsum)
    return S, freqs

def ft3(s, N):  # More efficient form of wikipedia equation

    S = []

    slice  = 0.0
    sliver = 2*np.pi/float(N)

    for k in range(N/2):

        sum_real = 0.0
        sum_imag = 0.0
        angle = 0.0
        for n in range(N):
            sum_real +=  s[n] * np.cos(angle)
            sum_imag += -s[n] * np.sin(angle)
            angle += slice

        slice += sliver
        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5
        S.append(tmpsum)

    return S

def ft4(s, N):  # Using wikipedia equation

    S = []

    for k in range(N/2):

        sum_real = 0.0
        sum_imag = 0.0
        for n in range(N):
            sum_real +=  s[n] * np.cos(2*np.pi*k*n/float(N))
            sum_imag += -s[n] * np.sin(2*np.pi*k*n/float(N))

        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5
        S.append(tmpsum)

    return S

def ft5(s, N):  # Roots of Unity Weighted Sum

    sliver = 2 * np.pi / float( N )

    root_real = np.zeros( N )
    root_imag = np.zeros( N )

    angle = 0.0
    for r in range(N):
        root_real[r] =  np.cos( angle )
        root_imag[r] = -np.sin( angle )
        angle += sliver

    S = []

    for k in range( N/2 ):

        sum_real = 0.0
        sum_imag = 0.0
        r = 0

        for n in range( N ):
            sum_real += s[n] * root_real[r]
            sum_imag += s[n] * root_imag[r]
            r += k
            if r >= N : r -= N

        tmpsum = np.sqrt( sum_real*sum_real + sum_imag*sum_imag )
        S.append( tmpsum )

    return S

def main():

    N  = 200
    fs = 100.0

    time_step = 1.0 / fs
    t = np.arange(0, N * time_step, time_step)

    f = 5.0
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs)

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S = ft3(y, N)
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    print len(S), len(freqs)
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()


RLS Algorithm (Memoryless) I have been studying the adaptive filters lately and now when I am at RLS (Recursive Least Squar) Algorithm I came across the term used in the weighting function of the RLS called forgetting factor ($\lambda$). This term ‘memoryless’ itself confuses me. My doubt is:

How is memoryless RLS different than the standard RLS? Does it involve some specific value for lambda (as in, it is said that if $\lambda$ $= 1$, system has infinite memory).

I would be glad if anyone could help.

In its bare classical form the RLS algorithm, recursively (for every new iteration), solves the classical problem of least squares; by computing the optimal FIR transversal filter coefficients $w[n]$ which minimizes the sum of error squares between the produced filter output $y[n]$ and a desired response $d[n]$ where the error sum is based on the entire data set beginning from the first sample $y[0]$ to the last sample (the current sample) $y[n]$. Therefore as time progresses, $n \to \infty$, the range of data set also goes to infinity, which means that the memory of the algorithm is infinite; i.e., the current value of the output $y[n]$ and the FIR coefficients $w[n]$ depend (indirectly via recursion) on the entire sample set.

A consequence of the infinite memory is that the estmation becomes more and more exact and refined given that the inputs were WSS random process. On the other hand a clear drawback of this infinite memory dependency is the fact that if the inputs are not WSS then the filter response to reflect the new conditions is slow (the entire past which is irrelevant to the current computation still effects the output due to infinite memory) and tracking is lost.

In order to improve this last situation, and therefore to make the algorithm track non stationary inputs, its memory should be shortened; i.e., only a finite amount of last local samples should effect the current computations. This is achieved by a forgetting factor $\lambda$ embedded into the recursion equations. The range of $\lambda$ is between $0$ and $1$ where $\lambda = 1$ converts the algorithm to classical form with infinite memory with no tracking and $\lambda \to 0$ has less and less memory and more responsive with better tracking...


What is the benefit of using symmetric kernel in Fourier transform? In image processing, a forward transform of an $M\times N\text{-}$pixel image $f$ from spatial domain coordinates $(x, y)$ to transform domain coordinates $(u, v)$ can be defined as [1, p.12]:

$$T(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)\cdot r(x,y,u,v)$$

The transformation kernel $r(x, y, u, v)$ is said to be symmetric if [1, p.13],

$$r(x, y, u, v) = r_1(x, u) \cdot r_1(y, v)$$

The forward kernel in discrete Fourier transform is [1, p.35]:

$$r(x,y,u,v) = e^{-j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)}\\(u=0,1,\dots,M-1,\,v=0,1,\dots,N-1)$$

My question is, What is the benefit of using symmetric kernels in Fourier transform?

I know the benefits of separable kernels, already.

[1]: "Fourier Transform" Slide set, George Bebis (UNR). Based on Ch 4 of “Digital Image Processing”, Gonzales and Woods. http://www.cs.umb.edu/~duc/cs447_647/spring13/slides/FourierTransform.pdf (archived copy)

A set of symmetric kernels along Cartesian axes should make the transform consitent with respect to (at least 90 degree) rotations of the input.  If the 2D input yields a particular 2D output, rotating that same 2D input by an angle of 90 degrees, should yield the same 2D output rotated by 90 degrees.  If the kernels in the two dimensions were not symmetric, the transform would yield different output that could not be rotated back to match the original output.

I'd have to think if the invariance applied to angles that were not multiples of 90 degrees.

You may also want to have a look at the Abel transform which has a circularly symmetric kernel, useful for physical situations with circular symmetry.


What is the meaning of imaginary values in the time domain? Say I take a real valued function in the time domain, take the Fourier transform to get its frequency domain representation, do some work on it in the frequency domain, which leads to imaginary values, and then take the inverse Fourier transform to get a new time domain signal.  What does it mean if the values of the time domain signal are imaginary?  I understand that imaginary values in the frequency domain allow for the representation of phase of each frequency (https://stackoverflow.com/questions/10304532/why-does-fft-produce-complex-numbers-instead-of-real-numbers).  But what do imaginary numbers mean in the time domain?

I can tell you what my particular application is, but I believe that this is a general question that others might come across in the future.

Imaginary values in the time domain simply mean that you actually have two signals, the real-part and the imaginary part. If you start out with a real-valued signal, take the Fourier transform, do some processing in the frequency domain, then transform back, and you end up with a complex-valued signal, then 9 out of 10 times you've done something wrong (or, if the imaginary part is very small, it is numerical noise due to round-off errors). In other words, in most cases the frequency domain processing should not make your originally real-valued time domain signal complex-valued, i.e. it should retain the conjugate symmetry of the signal's Fourier transform.


How can I take a fixed number of bins after N-point DFT when N is unknown? I am working with machine learning for time series classification. I am trying to extract features from the amplitude spectrum.

My current concern is that I cannot tell the length of the signal in advance. Here's one example signal with $N=192$ observations, corresponding to approximately 2 seconds of data collection:

I do not have an upper bound on the length of the signal, although more than 1,000 observations usually indicate I have an outlier.

My approach has been to take the input signal $x=(x_1,x_2,\ldots,x_N)$ and then find the one-sided amplitude spectrum from its N-point DFT:

$\mathcal{F}(x) = (X_1, X_2, \ldots, X_N)$
$A(X)=(|X_1|, |X_2|, \ldots, |X_{\frac{N}{2}+1}|$)

Given that my sampling rate is approximately $f_s=100\,\text{Hz}$, I figured I would interpolate the one-sided amplitude spectrum to 51 points. (PS: I made a mistake here; please see clarification below.)

This procedure gives me a 192-point DFT and a one-sided amplitude spectrum with 97 coefficients, which I interpolate to 51 (linear interpolation). For comparison, here's the original amplitude spectrum and its interpolation. The horizontal axis is given in Hz.

They don't look very different. However, the features I take from the spectrum are such as short-term energy, variance, standard deviation, number of peaks, and location of the first peak, as well as the magnitude of specific coefficients (am I right in saying that's the energy at those bins?). So different approaches in finding those 51 coefficients can lead to quite different results.

Also, I have the impression that interpolating the amplitude spectrum is just plain wrong. Looking at the figure, it seems obvious that the interpolation made me lose a lot of energy content in the range near 0-5 Hz.

For disclosure, I did check each related question suggested by Stack Overflow as I wrote mine. The most similar question to mine seems to be "How to combine bins of my DFT", however I believe I have several additional doubts---to start with, I'm eager for criticism on whether I am doing the right approach in taking the absolute value of the Fourier coefficients and then trying to reduce the number of bins, or whether I should zero-pad my signal to a number of observations that is a multiple of 50 (and thus hopefully avoid leakage? I don't know, I'm really struggling to understand the subject and I might be working with several misconceptions here; any corrections will be appreaciated)

PS: I made a mistake in my assumptions. I incorrectly assumed Matlab's interp1 function interpolates over a set of query points and averages the values in-between. What it is in fact doing is downsampling.

PPS: here's the code for coefficient reduction:

qp = (0:50) * (numel(X) - 1) / (50) + 1;
X50 = interp1(X, qp);

Selection of feature vectors is a bit of an art when a sufficient statistic isn't obvious.  If you want bins corresponding to the bandwidth of an equivalent fixed DFT.

I suggest you perform the STFT over your time duration varying signal, and then average the magnitude or if you want to avoid square roots the magnitude squared, over successive frames.  It will preserve "energy". This does complicate things by adding issues such as windowing and overlapping but sometimes having more parameters to tweak provides flexibility.

One way of looking at feature set selection is:

P. M. Baggenstoss, "Class-specific feature sets in classification," in IEEE Transactions on Signal Processing, vol. 47, no. 12, pp. 3428-3432, Dec 1999.
doi: 10.1109/78.806092
keywords: {probability;signal classification;statistical analysis;Class-specific feature sets;dimensionality problem;likelihood ratios;probabilistic classifiers;signal classification;sufficient statistics;Bayesian methods;Signal analysis;State estimation;Statistics;Training data;White noise},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=806092&isnumber=17473


Filtering without using any transforms A composite signal is defined as

$$x(t) = \sin(2\pi 10t) + \sin(2\pi 40t) + \sin(2\pi 60t) + \sin(2\pi 50t)$$

Separate the four frequencies without using any of the transform techniques.

You can achieve this result by using two combs filters : https://en.wikipedia.org/wiki/Comb_filter

Put simply, the comb filter consists of adding a delayed version of the signal to itself, causing destructive or constructive interference.

For instance, with $K = 20$ and a negative gain value after the delay line, you can significantly decrease or suppress harmonics that are even multiples of 20: $\sin(2\pi40t)$ and $\sin(2\pi60t)$, which leaves you with only
$ y(t) = \sin(2\pi10t) + \sin(2\pi50t)$.

Apply the same method one more time with $K = 100$ and you will end up with $y(t)=\sin(2\pi10t)$, although scaled and delayed due to the filters. The same method works for each of the four harmonics.

It might not be the most straightforward way to achieve this result, but it only uses simple filter structures. Another approach could be to use notch filters to remove the unwanted harmonics from your signal.


How can I change the formula to detect the barcode in every angle? - MATLAB How can I change the formula to detect the barcode in every angle?

 rgb = imread('barcode4.jpg');
% Resize Image
rgb = imresize(rgb,0.33);
figure(),imshow(rgb);
% Convert from RGB to Gray
Igray = double(rgb2gray(rgb));
% Calculate the Gradients
[dIx, dIy] = gradient(Igray);
B = abs(dIx) - abs(dIy);
% Low-Pass Filtering
H = fspecial('gaussian', 20, 10);
C = imfilter(B, H);
C = imclearborder(C);
figure(),imagesc(C);colorbar

If the filter you have works well enough for you, you could simply use imrotate to create several rotated versions of the source image, and use the filter you have on each of them (make sure to use bilinear to get sensible derivatives)


Transposed Direct Form II In the book digital Signal Processing by Proakis 4th Edition, at page 593 a table is drawn about some second order modules for discrete time systems. In that table there Transposed Form II drawn But when I searched the net, it was drawn completely opposite to that drawn in the book. On the internet, the input and output were placed on right and left respectively. I want to ask that which form is correct?

Direct-Form II is the canonical structure to implement an LCCDE, using minimum number of delay elements and multiplications. The input $x[n]$ is placed on the left, and the output $y[n]$ is on the right, complying with the signal processing conventions.
DF-II structure, implements the poles (IIR part) before the zeros (FIR part), however this may be undesirable for numerical issues with high gain poles.
The transposed DF-II structure, which maintains the canonical advantage of DF-II while avoiding the potential numerical instability of it, is obtained by reversing all branch directions, and replacing sum nodes with signal nodes, also replacing (and renaming) the input output nodes too; hence input becomes at the right side and output emerges at the left node. Transposed DF-II implements zeros first, and poles later.
However, for convenience, to comply with the convention that the input is applied at the left node, the transposed DF-II structure is horizontally flipped. This way, input $x[n]$ appears at the left again.
Both forms of transposed DF-II are equivalent.


How to compare two speech recognition softwares? I am working on building speech recognition software for my company.

In the meantime however, we will purchase some off-the-shelf commercial speech recognizers until the in-house one is ready. My question:
1. How are speech recognizers evaluated?
2. Are there any standardized tests available that can be used to test speech recognizers and determine the better one of the two?

NIST has a web page that provides a good detailed discussion of how to score ASR software.

http://www.isip.piconepress.com/projects/speech/software/tutorials/production/fundamentals/v1.0/section_04/s04_03_p01.html


Document image vanishing point detection and evaluation It is a natural idea to substitute our scanner with mobile phone, and perform OCR with the phone at the same time. If this can come true, two problems must be tackled: one is related to document image binarization, and the other is related to geometric distortion correction. Among all the distortions, projective distortion is most common, and in order to correct this kind of distortion usually the horizontal vanishing point and vertical vanishing point must be identified. There are many different ways to obtain the vanishing points, and different methods can lead to different vanishing points. Then my question is given several candidate vanishing points (horizontal vanishing points, for example) how can we select the best one? Of course we can get the projective transform model based on each vanishing point, and then correct the distorted image. Visual evaluation the corrected image can tell which vanishing point is better. But could I find the best vanishing point just before the correction? Thanks!

why use OCR techniques that are sensitive to vanishing points?

projective distortion just skews the letters- right?  So, why not just model those expected skews into the OCR templates for each letter?  Also, the nearby letters would skew similarly, so you could probably use correlated parallel lines among neighboring letters to estimate the amount of local/global skew (and avoid dealing with VPs).

also, I assume you mean the phone will take a pic of the document, so don't you also have to worry about the paper not being flat and having curve linear distortions too, which seems to me would completely confound, and supersede any vanishing point issue you are worried about.


How the FFT takes a cosine or sine and outputs the frequencies of the complex form? If i take the Fast Fourier Transform (FFT) of a cosine function, what has turned this cosine function into its complex exponential form which consists of $e^{i \omega t} + e^{-i \omega t}$ ?

Because on its own, cosine has a single frequency at omega, but the FFT must be interpreting cosine as a sum of its complex exponential because its output is amplitudes at frequencies $-\omega$ and $\omega$. So whats going on here?

I will try to explain it in an intuitive, not rigorous way.

The main idea behind the Fourier Transform is to "project" one signal $s(t)$ on another $\phi(t)$. If the projection is not zero, then $\phi(t)$ is "included" in $s(t)$, in the sense that you can write $s(t)$ as $$s(t) = A\phi(t) + \text{other signals},$$ with $A \neq 0$.

The projection is defined as $$A = \int_{-\infty}^\infty s(t) \phi^*(t)\,dt.$$

You can of course project $s(t)$ on multiple signals; you can find the projection on $\theta(t)$ $$B = \int_{-\infty}^\infty s(t) \theta^*(t)\,dt$$ and then you can say that $$s(t) = B\theta(t) + \text{other signals}.$$

In the case when the projection of $\phi(t)$ on $\theta(t)$ is zero (in this case they are said to be orthogonal), you can go further and say that
$$s(t) = A\phi(t) + B\theta(t) + \text{other signals}.$$ Since you found two components of $s(t)$, you can expect the "other signals" portion of this last equation to be "smaller" (have less energy) than when you project on $\phi(t)$ or $\theta(t)$ alone.

Now, what Fourier says is that if you project your signal on the set $\lbrace e^{j\omega t} \rbrace$, the "other signals" part of the equations above will be zero. So, this set is special: it's said to be a "complete" basis for all signals (it's not the only complete basis, though).

So, if you look at the Fourier Transform integral, you'll notice that it is projecting your time-domain signal on this complete set. It turns out that, when you project the cosine $\cos(\omega_0 t)$ on this set of exponentials, you get exactly two projections that are not zero. This means that the cosine can be written as the sum of these two projections.


What papers or resources explain the conditions needed for filter stability in real-time tuning of parameters? What papers or resources explain the conditions needed for filter stability in real-time tuning of parameters?
That is, how can I know when a filter is such or how to design a filter that will not pop or "jump discontinue" when its parameters are tuned?

Instead of trying to fit everything in a comment, I will answer the question here.
Assessing stability varies to which method you use to describe a filter (DE, Laplace domain, State - Space, FIR, etc). Stability itself is a very elaborate and interesting topic in controller and filter design, but I must keep it short enough. In essence, if you apply an impulse response to your filter, the energy of the output must either decrease or stay constant to make the filter (marginally) stable. But this does not tell you how this is linked to the parameters. Depending on your filters description, the following things should be checked for stability:

Laplace domain: roots of the denominator $\leq0$
Discrete Z-domain: absolute value of the roots of the denominator $\leq1$
Continuous State Space: Eigenvalues of A $\leq0$
Discrete State Space: absolute value of the Eigenvalues of A $\leq1$

So, suppose you have a simple LPF in the Laplace domain: $$F(s) = \frac{\omega_0}{s+\omega_0}$$
Then it can be seen that it is stable for any value of $\omega_0\geq0$
Next, lets take a standard FIR system in Z-Domain:
$$H(z) = \frac{b_0 + b_1z^{-1} + ... + b_Nz^{-N}}{1} = \frac{b_0z^{N} + b_1z^{N-1} + ... + b_N}{z^{N}}$$
Since the parameters do not affect the denominator, a FIR filter is always stable!

Is stable to parameter tunings the same as "being continuos function"

If your filter is linear, any continuous input will yield a continuous output, regardless of the parameters used. Any discontinuous input might yield an discontinuous output, but it remains stable. If your filter is non-linear, the methods to assess stability are much harder (Lyapunov equations is one I know), especially to link them directly to parameter tuning (you could trigger a bifurcation point where the filter is stable for some values, but unstable for others).
I hope this clarifies a bit and I am not currently explaining the obvious.


Bandwidth of a communication system As a rule of thumb, the bandwidth of a communication system with a complex signal must be equal to the bit-rate or symbol-rate?

Besides, after reading this post, it is not clear for me how a band-filter or low-filter signal can influence the maximum data rate.

In order to avoid inter symbol interference (ISI) the baseband transmit signal must meet Nyquist's first criterion. In frequency domain it can be formulated as follows.
$$
\omega_\mathrm{N} = \frac{v_\mathrm{s}}{2}
$$
With $v_\mathrm{s}$ the symbol rate and $\omega_\mathrm{N}$ the 3dB frequency of a low pass filter with symmetric flank. The flank must be point symmetric with regard to $\omega_\mathrm{N}$ and zero for $\omega > \omega_\mathrm{N}(1 + \alpha)$. Consequently, the baseband bandwidth $B=\omega_\mathrm{N}(1 + \alpha)$.

After modulation the bandwidth is doubled and the radio frequency bandwidth is $2B = 2\omega_\mathrm{N}(1 + \alpha) = v_\mathrm{s}(1 + \alpha)$. For the ideal case of a rectangular low pass filter $\alpha=0$ and $2B = v_\mathrm{s}$ which corresponds to the statement that the required bandwidth is equal to the symbol rate *). Similarly, if the bandwidth $B$ is given, the symbol rate is limited by $v_\mathrm{s} \leq \frac{2B}{1+\alpha}$

The required bandwidth indirectly depends on the bitrate, because a certain number of bits is transmitted per symbol.

The lowpass mentioned above is also referred to as impulse shaper and $\alpha$ as roll-off factor.

*) As an ideal lowpass isn't feasible, the required RF bandwidth in practice is always larger than the symbol rate.


Differences in transforming a sinc FIR filter from lowpass to highpass using spectral inversion and spectral reversal I have designed a lowpass sinc 1023 tap FIR filter with a Blackman window and this works fine.
I understand there are then two ways to transform this lowpass filter to a highpass one.
One is spectral inversion where you invert every coefficient and add 1 to the centre coefficient.
https://tomroelandts.com/articles/how-to-create-a-simple-high-pass-filter
The other is spectral reversal where you invert every other coefficient but only after generating a new lowpass filter with the frequency reflected about Fsample/2.
https://tomroelandts.com/articles/spectral-reversal-to-create-a-high-pass-filter
Both methods work well, giving similar though not quite identical results in the passband.  But the stopbands are quite different, with spectal reversal having much deeper rejection in the stopband.   However I cannot find any literature explaining why there is a difference.  Any suggestions on why the two methods are different or where I can read up more about them ?   Thanks in advance.

This answer is gonna come in installments.
First consider the ideal $\operatorname{sinc}(\cdot)$ filter, which in the frequency domain is an ideal brickwall.  We'll do this in normalized angular frequency, $\omega$.
Let
$$ H_\mathrm{LPF}(e^{j\omega}) = \begin{cases}
1, \qquad  |\omega| < \omega_0 \\
0, \qquad  \omega_0 <  |\omega| \le \pi \\
\end{cases} $$
and define the frequency response outside of $\pm\pi$ to be periodic:
$$ H_\mathrm{LPF}(e^{j(\omega+2\pi)}) = H_\mathrm{LPF}(e^{j\omega}) \qquad \forall \omega \in \mathbb{R} $$
This will have impulse response of
$$ h_\mathrm{LPF}[n] = \frac{\omega_0}{\pi} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) $$
Now, if you make a high-pass filter by subtracting the low-pass from a wire:
$$ H_\mathrm{HPF}(e^{j\omega}) = 1 - H_\mathrm{LPF}(e^{j\omega}) $$
Then the cutoff frequency remains the same $\omega_0$ and the impulse response is
$$\begin{align}
 h_\mathrm{HPF}[n] &= \delta[n] - h_\mathrm{LPF}[n] \\
 & = \delta[n] - \frac{\omega_0}{2} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) \\
\end{align} $$
Now, suppose you make a high-pass filter by heterodyning the low-pass from DC to Nyquist:
$$ H_\mathrm{HPF}(e^{j\omega}) = H_\mathrm{LPF}(e^{j(\pi+\omega)}) $$
Then the cutoff frequency becomes $\pi-\omega_0$ and the impulse response is
$$\begin{align}
 h_\mathrm{HPF}[n] &= (-1)^n h_\mathrm{LPF}[n] \\
 & = (-1)^n \frac{\omega_0}{\pi} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) \\
\end{align} $$
So, if we harmonize the two definitions of the cutoff frequencies, the two brick-wall high-pass filters will be identical in the frequency domain, so they should in the time domain.  But that would mean that
$$\begin{align}
\delta[n] - \frac{\omega_0}{\pi} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) &\stackrel{?}{=} (-1)^n \frac{\pi - \omega_0}{\pi} \operatorname{sinc}\left( \frac{\pi-\omega_0}{\pi} n \right) \\
\\
&= (-1)^n \frac{\pi - \omega_0}{\pi} \frac{\sin\left( \pi \frac{\pi-\omega_0}{\pi} n \right)}{\pi \frac{\pi-\omega_0}{\pi} n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{\sin\big( (\pi - \omega_0) n \big)}{\pi n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{\sin(\pi n)\cos(\omega_0 n) - \cos(\pi n)\sin(\omega_0 n)}{\pi n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{0\cos(\omega_0 n) - (-1)^n\sin(\omega_0 n)}{\pi n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{- (-1)^n \ \sin\left(\pi\frac{\omega_0}{\pi} n\right)}{\pi n} \qquad & (n \ne 0) \\
\\
&= \frac{- \sin\left(\pi\frac{\omega_0}{\pi} n\right)}{\pi n} \qquad & (n \ne 0) \\
\\
&= -\frac{\omega_0}{\pi} \frac{\sin\left(\pi\frac{\omega_0}{\pi} n\right)}{\pi \frac{\omega_0}{\pi} n} \qquad & (n \ne 0) \\
\\
&= -\frac{\omega_0}{\pi} \operatorname{sinc}\left(\frac{\omega_0}{\pi} n \right) \qquad & (n \ne 0) \\
\end{align}$$
So I guess it's true, because you can also show that equally exists for the case of $n=0$.


python scipy fft on numpy hanning window smears peaks I have a signal that I have acquired from an experimental instrument, that I wish to examine in the frequency domain.  I don't care about phase for this exercise, I only care about magnitude.  I deliberately planned my experiment so that the dominant frequency would fall completely within a single fft bin.  when I use the scipy fft function on an unfiltered window, the fft shows a clean spike as expected.  However, when I first apply a numpy.hanning window, the spikes become smeared.
Note the mean of the signal (the zero bin) also shows the same smearing effect.
My code is as follows:

eta=[instrument data series of length 2048]
etaHann=np.hanning(nfft)*eta
EtaSpectrum=abs(sp.fft(eta))
EtaSpectrum=EtaSpectrum*2/nfft # convert to amplitude
EtaSpectrumHann=abs(sp.fft(etaHann))
EtaSpectrumHann=EtaSpectrumHann*2*2/nfft # also correct for Hann filter
frequencies=np.linspace(0,samplingFrequency,nfft)
plt.plot(frequencies[0:200],EtaSpectrumHann[0:200],'b.-',label='Hann filtered')
plt.plot(frequencies[0:200],EtaSpectrum[0:200],'c.-',label='rectangular')
plt.xlim(-0.01,2)
plt.xlabel('frequency')
plt.ylabel('amplitude')
plt.grid()
plt.legend(loc=0)
plt.show()

Plotted outputs below:

Is this meant to happen?  Does anyone else get this problem?  For my data, it seems a rectangular window is fine, but to me it seems that the Hanning window must be wrong here, as the integral of the entire spectrum should be the same (after I multiplied by 2 to correct amplitude), but in this case it seems to be different.

Is this meant to happen?

Yes, that is absolutely a well known effect of using any window function.
Taking a look through the Wikipedia article on window functions, we find that the rectangular window has the sharpest peak.
This is quite simply due to the fact that the rectangular window acquires the most data, meaning that it can distinguish between slightly different frequencies the best.

The problem with the rectangular window is that it has large side lobes.
In other words, if you measure a signal which doesn't exactly match one of the discrete Fourier transform (DFT) frequencies (a.k.a. bins), then with the rectangular window you get a lot of leakage into neighboring bins.
The formula for this is actually rather simple.
Suppose your sampling interval is $\delta t$ and you measure $N$ points.
Then the DFT frequencies are
$$f_n = n / N \delta t, \quad n\in[0,1,\ldots,N-1] \, .$$

Now suppose you measure a sinusoid with frequency $\xi$ so that the sampled values are $x_n = \cos(2 \pi \xi n \delta t)$.
The mod squared of the DFT values are then
$$ |X_k|^2 = \left( \frac{\sin(\pi(\xi - k))}{\sin(\pi(\xi-k)/N} \right)^2 \, . \tag{$*$}$$

We plot this function in Figure 1 for $\xi=24.4$.
Note that when $\xi$ is not an integer the DFT values are nonzero for all bins, and you can see from equation $(*)$ that the values of $|X_k|^2$ only fall off as one over the square of frequency.

Figure 1: DFT values for a cosine of frequency 24.4.
The blue line indicates the underlying function given in equation $(*)$ and the black dots are the values of that function evaluated at the integers.
In other words, the black dots are the DFT bin values.

The point of a window function is to make the leakage of power into bins near the real frequency (in this case 24.4) fall off faster, but this comes at the cost that the main peak is broadened.$^{[a]}$

  Does anyone else get this problem?

It's not a problem.
It's how windows are supposed to work.

  For my data, it seems a rectangular window is fine, but to me it seems that the Hanning window must be wrong here, as the integral of the entire spectrum should be the same (after I multiplied by 2 to correct amplitude), but in this case it seems to be different.

Applying a window changes the total power in your signal, as you have observed.
This is not a mistake at all.
It's pretty obvious that this has to happen because the integral (sum) of the power in the frequency domain has to equal the power you would have computed in the time doman (Parseval's theorem).
It's clear that multiplying e.g. a cosine by a window function changes the integral of the time domain signal.
This is called "loss" and must be taken into account when using windows!

$[a]$: I'm not actually showing this because I didn't bother to compute the points with the window, but the Wikipedia article shows many examples.


What is the meaning of the length of diagonal lines in FFT I am relatively new to FFT analysis of images and have a question about them. I like to think I understand the basics but many of the properties shown in a FFT image are beyond me.

My main question is what is what is the meaning of the length of diagonal lines in a fourier transform? What is the meaning of the length of a diagonal line on the image domain or in the Fourier domain? The diagonal lines are related to the edges in the image but what kind of property does the length describe?

If asking these questions means that I don't fully understand FFT, are there any resources that would help me dig deeper?

If you rotate an image (2D function), its 2D Fourier transform also rotates.

Those lines mean you have something like a barcode oriented in that direction in your image. In one direction it has lots of variations and behaves like a white noise with flat spectrum and in the other direction it has very slow variations.

I'm not sure what do you mean by the length of diagonals, but it's a visual measure of maximum frequency available in those directions.


Examples of mixture of 'almost periodic' signals/data in nature or practice I am working on algorithms to separate mixtures of 'almost periodic' signals. One such example is composite ECG consisting of maternal and fetal signal.

Can anyone help with examples of such signals? The signal can also be from internet data.

The sound produced using the lowest pitched strings of large stringed music instruments often includes slightly inharmonic overtones (due to finite string thickness and stiffness), which renders the waveform of the note pitch not quite, but "almost periodic".  You can record these sounds from the bottom octaves of most (physical, stringed) pianos, and bass guitars, etc.


What can we understand when we increase the $N$-point DFT for a given signal? I have an input sequence $$x(n)=\cos(0.48\pi n)+\cos(0.52\pi n).$$ I am determining its spectrum (amplitude of DFT values) based on the finite number of samples.

For Example :

Taking the first $10$ samples of $x(n)$, $n=0, \ldots, 9$, calculating the $10$-point DFT of $x(n)$
Using the $10$ samples from part $1$, calculating and plotting the $100$-point DFT.
Taking the first $100$ samples of $x(n)$, $n=0, \ldots, 99$, ploting the $100$-point DFT of $x(n)$.

Questions:

Why are we going for a $100$-point DFT/ $200$-oint DFT ?
And what happens when we take $10$ samples and calculate $100$-point DFT ?

If the signal length $M$ remains the same but only number of points $N$ in DFT is increased (assuming $M \leq N$) then this will give you the finely spaced samples $X[k]$ of the same true spectrum $X(e^{j\omega})$ as given by the DFT frequency sampling relation: also known as the zero-padding result
$$X[k] = X(e^{j\frac{2 \pi}{N}k}) = \sum_{n=0}^{N-1} {x[n]e^{-j\frac{2 \pi}{N}k}} = \sum_{n=0}^{M-1} {x[n]e^{-j\frac{2 \pi}{N}k}} ~~~ , ~~~~  k=0,1,2,...,N-1 $$

On the other hand, if both the analysed signal length $M$ and the DFT length $N$ are increased (assuming $M=N$) by using a longer analysis window (considering $x[n]$ defined for a very long, possibly infinite length duration) then the windowed true spectrum of the underlying signal whose samples are investigated, is also enhanced, thereby providing samples $X[k]$ of a higher frequency resolution spectrum of the very long length signal $x[n]$


Is it necessary to downsample data for filtering of low frequencies? My data is sampled at $f_s = 32kHz$, and for a particular analysis I'm only interested in signals with certain frequencies from $f_{highpass} = \text{0.1Hz}$  to $f_{lowpass} = \text{300 Hz}$.

When I build and apply a 4th order butterworth filter to my data I get an error and there's not much left of the signal. Here's some Matlab code:

order = 4;
[b,a] = butter(order, [0.1 300]/(32000/2), 'bandpass');
y = filtfilt(b,a,data);

The error message I get is: Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.101227e-16.

I've also tried to apply a low pass filter first, then followed by a highpass filter:

FOI = [0.1 300]
SampleRate = 32000;
[b_l,a_l] = butter(order, FOI(2)/(SampleRate/2), 'low' );
[b_h,a_h] = butter(order, FOI(1)/(SampleRate/2), 'high' );
tempData = filtfilt(b_l, a_l, data);
output = filtfilt(b_h, a_h, tempData);

But I keep getting the same result.

My questions:

why does the filtering fail?
why does the filtering with order 3 work but not with order 4?
do I need to downsample the data (e.g. to 1kHz) before I can do the filtering?

EDIT: other people have hit the same problem but I'm looking for a definite answer.

EDIT2: if I change the highpass value to 1 and do the highpass filtering before the lowpass filtering, everything works. But why? Why is the cutoff of 0.1 causing problems?

Yes, you can get much better results if you downsample your data first. The problem occurs because the width of your transition band relative to the original data sample rate is very small. This can cause numerical problems when designing or implementing digital filters. The problem only gets worse as the filter order increases, especially with IIR filters.

With that said, in some cases, you might have more success (i.e. better numerical stability) if you convert the filter to second-order sections. However, for this case, I think the right answer is to downsample the data first.

The problem likely occurs because of the very small cutoff frequency on the low end of your bandpass filter. There is very little room for the transition band (less than 0.1 Hz, which is 1/320000-th of the sample rate), which is difficult to implement. If you downsample the data as much as possible, then you can at least increase this ratio by a good amount, which should give you better results.

I would also try using FIR filter implementations if possible, as they typically have superior numerical stability due to their lack of feedback.


FFT spectrogram in log frequency space. How? I was reading a paper and the authors say that they extract FFT spectrograms and log-scale in frequency domain. I understand how the general spectrogram is computed, but how to make it log-scale in frequency domain? Or maybe there are libraries in Python, which already implement this one?

You can use either interpolation between FFT bins, or something like an MFCC triangular filter bank, except using equal delta_log spaced triangles instead of Mel-spaced triangles, or both.


Beam resolution improvement of linear array using MVDR beamformer One of the major advantages of using MVDR beam former is to improve bearing resolution as compared to conventional delay & sum beam forming. The bearing resolution is defined by 3 dB beam width of linear array beam pattern. There are also formulas for calculating 3 dB beam width of arrays.

Is there any way to calculate the 3 dB beam width of linear array using ideal MVDR beam former (or beam pattern calculated by ideal MVDR) so that the bearing resolution improvement with respect to conventional beam former can be quantified for same length of array?

The MVDR beamformer is an adaptive technique - the response is dependent on the the interfering signals i.e. their locations, their SNRs and their correlations. This makes it difficult to say what the resolution is in traditional 3dB width type of terms. Essentially the MVDR will try to place nulls in the beampattern where the interferers are - so a 3dB width doesn't really tell you anything about the resolution.

You can also get strange effects, e.g. you'll get 0 dB gain in the desired look direction, and you'll have a sidelobe level of +10 dB somewhere, because it's trying to place a null in the beampattern in a particular location.

As a purely hypothetical example - consider you're forming a beam at 90$^\circ$ look direction and you have an interferer at 135$^\circ$. You apply the MVDR and you'll see your beampattern has a null at 135$^\circ$ and the mainlobe around 90$^\circ$ will have a certain width. Now if you add another interferer  at 45$^\circ$ and apply the MVDR - you'll have nulls in both the interfering bearing locations, but your mainlobe width around 90$^\circ$ has now increased. So that mainlobe width in the single interferer case doesn't really tell you much about what's going to happen when you add another interferer. Essentially you are using up another degree of freedom to null out the additional interferer.

What is probably more interesting to look at is - how close can I place an interferer before the MVDR can no longer null it out.

A good discussion on MVDR is in "Array Signal Processing" by Johnson and Dudgeon, Prentice Hall.

You can also look at "Optimum Array Processing" by Van Trees, Wiley Interscience - this may be overkill for what you want. You also have to be careful - he a uses slightly different definition of MVDR.


LTI system without constant coefficient differential equation I have encountered a system where the output $y(t)$ and input $x(t)$ are related in the laplace domain as:

$$Y(s) = H(s)X(s) \tag1$$

which is typical. However, $H(s)$ is not a rational function of polynomials of $s$. Instead, $H(s)$ contains terms like $\log(s)$.

As far I knew, one can write equation $(1)$ only for an LTI system and LTI systems have rational polynomial transfer functions. Am I correct in this analysis?

If I have a logarithm in $H(s)$, that implies that I cannot write a linear differential equation with constant coefficients in the time domain. Is it possible to have an LTI system that is not representable by such differential equation?

Does a frequency domain analysis always imply an LTI system?

As long as a system can be characterized by its transfer function, or, equivalently, by its impulse response, the system is LTI. For the LTI property it's irrelevant if the transfer function is rational or not.

If the transfer function is rational (and positive-real), the system can be realized by a circuit containing lumped elements, i.e., ideal R, L, and Cs. On the other hand, a simple example of a non-rational transfer function is an ideal delay line (e.g., approximated by a long cable), the transfer function of which is $H(s)=e^{-sT}$, where $T$ is the delay.

In sum, the exact expression for $H(s)$ is irrelevant for the LTI property of the corresponding system. However, $H(s)$ can only be realized by resistors, inductors, and capacitors if it is a rational function (and if it is positive-real).


Is a Fourier transform a sound way to analyse a transient signal? I am currently working on a project that involves analysing transient signals from sensors.

While not actually part of the analysis itself, I discussed it with the team, and they are using an fft to analyze the transient measurement.

From what i recall , Fourier Transform is a tool for periodic signals, not transient. Am I right? Am I missing something?

If the signal is zero-valued outside the time window being analyzed ("analysis window"), then the discrete Fourier transform (DFT) exactly samples the discrete-time Fourier transform (DTFT) of the signal at harmonic frequencies. DTFT implies no time-domain periodicity. The sampling of DTFT by DFT is evident from definitions of the two transforms, if we apply a finite support signal constraint:
$$\text{DTFT:}\quad X(\omega) = \sum_{n=-\infty}^{\infty} x[n] \,e^{-i \omega n}$$
If $x[n] = 0$ for all $n$ such that $n < 0$ or $n \ge N$ (finite support signal constraint), then:
$$= \sum_{n=0}^{N-1} x[n] \,e^{-i \omega n}$$
The frequency variables of the two transforms are related by $\omega = \frac{2\pi}{N}k$, so:
$$= \sum_{n=0}^{N-1} x[n] \,e^{-\frac {i 2\pi}{N}kn}$$
This is identical to the definition of DFT:
$$\text{DFT:}\quad X[k] = \sum_{n=0}^{N-1} x[n] e^{-\frac {i 2\pi}{N}kn}$$
Choosing a longer window (larger $N$) enables DFT to sample DTFT more densely.
For a signal that has non-zero values outside the analysis window, doing the analysis on just the window is equivalent to first multiplying the signal by a rectangular window, and then sampling the DTFT of this windowed signal using DFT. There are other window functions that may have more favorable characteristics for the type of signal and for the particular analysis task, allowing to isolate a transient of interest without too much distortion in time domain or smearing in frequency domain.
Fast Fourier Transform (FFT) is often used because it is very fast to compute, even though a decomposition into harmonic sinusoids might not always be the best way to describe a signal.


MATLAB: Why my Pixel Indices are so different from that in Spatial Coordinates The question is revised to make it more specific and constructive.

I am doing imtransform, yet the xdata and ydata are very much different from that of size(output).

[ output, xdata, ydata ] = imtransform( imgDouble{1}, tform );

Real data test:

K>> size(output)

ans =

         854        1230

K>> xdata(2) - xdata(1)

ans =

   1.6635e+03

K>> ydata(2) - ydata(1)

ans =

  876.1373

It is about Expressing Image Locations.

There are two kinds of image location expression,

One and the most common one is the

pixel indices:

in row and column

The other one is Intrinsic Coordinate System:

express in x and y:


Effect of cyclic prefix and zero-padding in OFDM What happens if zeros are added instead of cyclic prefix in OFDM system?
How does it affect performance of system?

The main idea of using Cyclic prefix in OFDM systems is because of 2 reasons,

1) To eliminate ISI from previous symbol.

2) To model linear convolution as circular convolution

If an OFDM symbol is Zero padded instead of cyclic prefix, ISI due to adjacent symbol still gets eliminated, but it cannot be modeled using  circular convolution. Circular convolution converts frequency selective multipath channel into flat fading channel. Circular convolution simplifies channel estimation and equalization which can be done easily in Frequency domain.


LTI system input upsampling Let's assume that a linear and time-invariant system is sampled at 2 different frequencies $F_{s}$ and $2F_{s}$ (e.g. 5Hz and 10 Hz). It gives
$$Y_{F_{s}}(z) = H_{F_{s}}(z)X_{F_{s}}(z)$$
$$Y_{2F_{s}}(z) = H_{2F_{s}}(z)X_{2F_{s}}(z)$$
where $Y_{F_{s}}$ and $Y_{2F_{s}}$ are respectively the z-transforms of the sampled signals, $H_{F_{s}}$ and $H_{2F_{s}}$ are the impulse responses of the systems and $X_{F_{s}}$ and $X_{2F_{s}}$ are its inputs. The signals $Y$ are known, the impulse response are also known :
$$ H_{F_{s}}(z) = \frac{a+bz^{-1}}{1+cz^{-1}+dz^{-2}}$$
$$ H_{2F_{s}}(z) = \frac{a'+b'z^{-1}}{1+c'z^{-1}+d'z^{-2}}$$
and the inputs are unknown.

I am wondering whether it is possible to express $X_{2F_{s}}$ in terms of $X_{F_{s}}$.

Since I am dealing with upsampling, I thought that $X_{2F_{s}}(z)=X_{F_{s}}(z^2)$ and so I had
$$Y_{2F_{s}}(z) = H_{2F_{s}}(z)X_{F_{s}}(z^2).$$ But I performed numerical simulations and it seemed to be false. So, where is my mistake and is it possible to express $X_{2F_{s}}$ in terms of $X_{F_{s}}$.

Any suggestions will be appreciated.

I am wondering whether it is possible to express $X_{2F_s}$ in terms of $X_{F_s}$.

In general, no.

This is the sampling theorem: unless you've bandlimited $X_{2F_s}$ first to be alias-free sampled by $F_s$, then the digital signal obtained by sampling twice as fast will contain information that the digital signal sampled at $F_s$ does not contain. And information that is not contained cannot be recovered.

If you, however, did not band limit $X$ to the Nyquist bandwidth dictated by the sampling rate, your $X_{F_s}$ will already contain aliases. Then you can simply look at the "upper" half of $X_{2F_s}$, shift it down to overlay exactly the lower half, which is exactly what's called aliasing, and reconstruct $Y_{F_s}$ from $Y_{2F_s}$.


How to extract the vein pattern from an IR image of finger? I am trying to isolate (or segment to be precise) the vein pattern from Infra red image of a finger (using MATLAB) like shown below:

So my output will be a binary image like (not the exact dimension of image shown):

I tried histogram equalization to make the veins more visible, but it was not that successful. Any help to proceed further would be very helpful.

You might want to also try running a ridge filter on the dataset; something like the frangi filter might work.

A sample result is

NB
the white border is due to the border in the image I pulled from this site


Visualizing the "frequency shift" theorem in FFT I was trying to visualize the frequency shift theorem in MATLAB, which states that if the Fourier transform G(f) is shifted by a constant $f_o$, its inverse transform  g(t) is multiplied by $e^{j2\pi tf_o}$.

Say one makes two symmetric real valued Gaussians (G(f)) in the frequency domain as shown in the code below as a function of bin numbers. The peaks are symmetrically displaced from the center bin by a constant. When the inverse FFT is applied G(f), the output g(t) looks look the attached photo. A displaced Gaussian in the frequency domain should produce a Gaussian multiplied by cosine waves.
Intuitively, it does appear like an exponential decay times a cosine wave.

a) The FFT of a Gaussian is a Gaussian, but g(t) looks like exponential decay times a cosine.
b) Another point is why do we see a two-sided symmetric output in the time domain? Shouldn't we see only one sided result in the time domain?

What is the discrepancy here in the time domain?

Thanks.

Frequency Shift Theorem

W=[1:1:9001]'; % bin numbers

n=1; % Gaussian shape, n=1 Gaussian

H=10; % Height of Gaussian

s=25; % standard deviation in frequency domain

w_left=1400; % Position of Gaussian as bin numbers

w_right=length(W)-w_left; % Position of Gaussian as bin numbers

G_left=H*exp(-((W-w_left).^2/(2*s^2)).^n);

G_right=H*exp(-((W-w_right).^2/(2*s^2)).^n);

GShifted=G_left+G_right;

Inverse_GShifted=real(ifft(GShifted));

figure (1)

subplot (1,2,1)
plot (W, GShifted,'r', 'LineWidth', 1.5);
title ('Shifted Gaussian in Frequency Domain')
xlabel('Bin Numbers (Frequency)')
ylabel('Real')
xlim([1 9001])

subplot (1,2,2)
plot (W, Inverse_GShifted,'b', 'LineWidth', 1.5);
title ('Inverse FFT of Shifted Gaussian')
xlabel('Bin Numbers (Time)')
ylabel('Real Part')
xlim([1 9001])

In the line

Inverse_GShifted=real(ifft(GShifted));

I have changed to

Inverse_GShifted=real(fftshift(ifft(GShifted)));

This is because the ifft output is from $0 \le n \le N-1$, but your actual time output is symmetric about $n=0$. So what you have to do is to center it around $n=0$ by doing fftshift() command. It will give the output of $-N/2 \le n \le N/2-1$ whose plot is below

Regarding the zoomed output, the IFFT of gaussian function $G(f)$ is again Gaussian $g(t)$. It is getting modulated or point-by-point multiplied by $e^{j\omega_0 n}$ where $\omega_0 = 1400/4500\pi$ (considering $-\pi \le \omega_0 \le \pi$), which is the shift in frequency that you applied. ie, $g[n].*e^{j\frac{14}{45}\pi n}$. The period of this waveform turns out to be $\approx 6.43$ so you will see the peaks repeated after every 6 or 7 samples which is correct as per the zoomed plot. $14/45$ being irrational may be the reason why period is not an integer.


SNR calculation in audio reconstruction and what is an acceptable value I have the original signal, a .wav file and by blind source separation, i have reconstructed the signal.
How can I calculate the SNR in this case? or is it good to use any other performance parameters? What should be an acceptable value of SNR?

The common practice in the source separation community is to report 3 metrics defined in Vincent et al's paper, called SDR (signal to distortion ratio), SIR (signal to interferences ratio), and SAR (signal to artifacts ratio).

The reason for this is that in source separation applications, there can be two types of "noise" in the reconstructed signal: noise due to mis-separation, called interferences (for example, if you want to extract the vocals from a song, this would be a residual of the background music), and noise due to the reconstruction algorithm itself, called artifacts (for example, glitches due to the STFT phase estimation process). The former is sometimes referred to as "crosstalk" in other contexts.

There is usually an application-dependent trade-off between the two - for example in a Karaoke application (vocal removal), it is not a big problem if the audio accompaniment has still a few traces of the original voice (low SIR), as long as the result is pleasing to the ears and glitch-free (high SAR). Most separation algorithms have "knobs" that will allow you to balance SAR and SIR. For example, increasing the STFT window size tends to increase SIR but decreases SAR.

The estimation of these metrics is described in the paper by Vincent et al, and implemented in the bss_eval toolbox. Of course, access to the original individual sources you want to recover is necessary! You will have to evaluate your algorithm on a data-set for which the "right answer" (original individual track) is known.

Regarding the acceptable value of the metrics, you need to compare your results with those published in the academia for the same problem you are attempting to solve. Keep in mind that these are objective metrics, and as such, do not take into account the perceptual quality of the noise. It might be possible that an algorithm with a -12dB SAR will be more pleasing to the ears than an algorithm with a -24dB SAR if the former produces artifacts well spread in time and below the masking threshold of the reconstructed signal. SDR/SIR/SAR can be used in development to rule out bad ideas and tune your algorithms, but for a commercial application, listening tests are ultimately necessary.


Understanding Parseval's Theorem with Discrete Wavelet Transform I have difficulty to understand the results I get with implementing Parseval's Theorem in Python to DWT. I have the good results getting the Energy with Fourier transform and the time series in python:
# Parseval theorem energy
def ParsevalTheorem(data):
    energy_sum = 0
    for i in range(len(data)):
        energy_sum += abs(data[i])**2
    return energy_sum

# dwt_data[0] => approximation component at final level, dwt_data[1:] => detail components
def DWTParseval(dwt_data):
    details_sum = 0
    for i in range(len(dwt_data)-1):
        details_sum += ParsevalTheorem(dwt_data[i+1])
    approx_sum = ParsevalTheorem(dwt_data[0])
    final_sum = approx_sum + details_sum
    return final_sum

fourierTransform = np.fft.fft(short_signal)
print("fourier energy: ", ParsevalTheorem(np.abs(fourierTransform))/len(fourierTransform))
print("Org energy: ", ParsevalTheorem(short_signal))
print("DWT energy: ", DWTParseval(app1)) # app1 is haar discrete wavelet transform using pywt.wavedec(data, "haar", level = 3)

Results:
fourier energy:  1305035.7546624008
Org energy:  1305035.7546624022
DWT energy:  1309077.6827128115

I've gathered the information on using Parseval Theorem from equation: Equation Link1

I have also encountered another equation to get the Energy but if I divide the Approximation sum with it's length it's in whole different scope than the original signal energy: Equation Link2

I some what understand the Parseval theorem when dealing with fourier transform, but lost with these equations when dealing with DWT.
PS: I know there is more Pythonic way to do the code but I intend to apply it in a different language also.

Parseval's identity and Plancherel's theorem finally boil down to orthogonality. When one decomposes a data (with samples), via a scalar product, onto an  orthogonal sequence (yielding coefficients), there exists a certain preservation (equality, up to a proportionality factor) of energy between samples and coefficients. There are some technical conditions, and in certain cases, one only get inequalities (cf. Bessel's inequality) or frame bounds.
The equation for the discrete wavelet transform (DWT) might be incomplete, with respect to the indices. I for instance think that in the second term of the RHS, the scaling factor should be $N_j$, not $N_J$ (and this depends a bit on how discrete wavelets are implemented). Basically, an orthogonal wavelet transform projects data onto basis elements gathered in groups called subbands. Each wavelet subband comes from $N_j$ vectors, with an additional $N_J$ vectors for the approximation. And normally, the total number of vectors  should be (about, honestly, this depends on signal extension) the number of samples $N$, in other words: $N=N_J +\sum_{j=1}^JN_j$.


What is the advantage of multiple antenna in receiver in energy harvesting ?harvested more power? What is the advantage of multiple antenna in receiver in energy harvesting ?

Assume there are $3$ transmitted antennas in the A device,and $1$ received antenna in the B device,now assume A device use beamforming with $3$ antenna to transmit some signal to B device,and B device harvest $0.1W$ power from the signal.

Now if i add one more received antenna to device B,that is, there are $2$ antennas in B device now,and assume these two antennas are very closed,how much power will device B harvest now? still $0.1W$? or $0.1W\times2$ ?

I am confused about this question,because if device B still harvest $0.1W$,then what is the advantage  of multiple antenna in receiver in energy harvesting ?

If device B harvest $0.1W\times2$,and this mean if i want to create power,i just use $1$ antenna to transmit 1W RF power to a receiver with $\infty$ antennas receiver,then i can get $0.1W\times \infty $ power ,but it is impossible .

After thinking,i think when the antenna number of device B become 2,the harvested power of B will not become $0.1W \times 2$,WHY?

Because when there is just one antenna in B device,the beamforming direction just aim that antenna,because just that antenna has to receive or harvest the signal,so that antenna can harvest all the RF energy

However,if the antenna number become two,device A has to let both of two antenna of B receive or harvest the signal,so the beamforming direction has to change,maybe the middle of these two antenna,so now first antenna will harvest about $0.05W$,and the other will receive $0.05W$ as well.Or perhaps $0.06W$ and $0.04W$,but the total power will not be larger than 0.1W


Optimal Level of Wavelet Decomposition for Denoising If a certain mother wavelet is defined and the optimal threshold shall be defined via Donoho, how can I define the optimal decomposition level to achieve the best denoising results.

I'd say it depends on the noise properties and of course the image itself.

What you can think is that most Denoise Filters can handle only the High Frequencies of the noise.
Hence the decomposition process moves Low Frequency of the noise to the High Frequency part for the spectrum.

So if your noise is white you need to go down as you can.
If it colored and most of its energy is in the High Frequency you can do only few levels.


why is the z transform transfer function 1/(z-1) called an integrator? I am reading up on delta sigma modulators and there this term $\frac{1}{z-1}$ that appears repeatedly and is referred to as an "integrator".  Why is this so ?

There are a couple reasons.  One is that $(1-z^{-1})$ represents $x[n]-x[n-1]$ which is a finite difference over a very small period of time.  and that is an approximation to a differentiator.  The reciprocal is

$$ \frac{1}{1-z^{-1}} = \frac{z}{z-1} $$

which is the inverse operator.  We normally call the inverse operation of differentiation, we call that "integration".

Another reason is simply to implement that term as a transfer function of a tiny little LTI system:

$$ \frac{Y(z)}{X(z)} = \frac{1}{z-1} = \frac{z^{-1}}{1-z^{-1}} $$

or

$$ Y(z)(1 - z^{-1}) = Y(z) - Y(z) z^{-1} = X(z) z^{-1} $$

that translates to

$$ y[n] - y[n-1] = x[n-1] $$

or

$$ y[n] = y[n-1] + x[n-1] $$

so the current output sample is the previous output added to the (slightly delayed) input.  the output is an accumulation of the input.  similarly an integral of an input $x(t)$ is an accumulation of that input until the present time $t$.

i would say a better representation of a discrete-time integrator is $\frac{z}{z-1}$.  that would correspond to

$$ y[n] = y[n-1] + x[n] $$


How to Increase the Resolution of a Video from a Sequence of Photos? I have a video filmed in a relatively low quality and resolution and a sequence of photos of higher quality taken of the same scene at the rate of about one image every two seconds or so. Could those photos be used to filter the video sequence to upscale it without loss to the sharpness?

One way to do it is to solve a MAP problem of the up scaled video and using the High Resolution images as a prior.

Try looking at the articles - Super Resolution MAP.


Remove audible click sound from high frequency hopped sound I have a total sound duration of $T$.

I break the duration to several sub-durations and in each sub-duration, I generate several discrete sinusoidal values by the formula $\sin(2\pi ft)$ with a sampling rate of $f_s$.

The frequency used in consecutive sub-durations is different.

Now when I generate sound using these values (using matlab) I get audible click sounds between the point where two sub-durations meet; I think because of discontinuity. How can I remove these audible click sounds?

Edit
the code I am using (roughly)

duration = 0.25;
sampleRate = 44100;
numberOfSamples = (duration * sampleRate);
frequency = 19000:500:21500;
freqHopSamples = [];

for j = 1:length(frequency)
    startIndex = ((j-1)*(numberSubSamples) + 1);
    endIndex = j*numberSubSamples;
    t = (startIndex:endIndex)/sampleRate;
    freqHopSamples(startIndex:endIndex) = sin(2*pi*frequency(j)*t);
end

wavwrite(freqHopSamples, sampleRate, 'sound.wav');

One method is to adjust the phase of the sinusoid at the beginning of each segment to minimized the 1st or 2nd delta phase difference from the phase at the end of the last segment.  Another method might adjust several samples at each boundary to minimize multiple derivatives.

You currently seem to start with a phase of zero, which could randomize the phase delta between segments, depending on frequency and segment length.


Gain error when performing STFT spectogram of increasing frequency signal I'm performing a STFT of a sine wave in which the frequency increases (up-chirp signal) using python numpy.fft as shown in the following code:

import numpy as np
import matplotlib.pyplot as plt

def stft(x, L, overlap, window):
    opoints = int(L*overlap)
    X = np.array([np.abs(np.fft.fft(window*x[i:i+L])[:L//2])
                 for i in range(0,len(x)-L,opoints)])
    return X

t = np.linspace(0, 5, 11025*5)

xf = np.linspace(0, 2000, 11025*5)
w = 2*np.pi*xf
x = np.sin(w*t)

f = np.linspace(0, 11025, 1024)

window = np.hanning(1024)

stft = stft(x,1024,0.5,window)

plt.imshow(stft[:,:1024//2].T, cmap = plt.get_cmap('gist_yarg'), origin='lower', aspect='auto',
                 interpolation='nearest', extent=(0,5, 0,11025/2))

plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.show()

It works for a constant frequency sinewave. However it shows not an offset but a gain error if using the up-chirp signal as input.

STFT of increasing frequency sine wave

I'm using:

Hanning window
Overlap = 0.5
Window size 1024
Sampling frequency 11025Hz
Total sample time: 5 seconds

What is causing the gain error that leads the 2 kHz signal appearing at double frequency at the end of the time window?

Thank you in advance,
Pedro

A sinusoidal signal is

$$y(t) = sin(\phi(t))$$

The phase of a linear chirp, is

$$\phi(t) = \phi_0 + 2\pi\left(f_0t+\frac{k}{2}t^2\right)$$
where
$$ k = \frac{f_1-f_0}{t_1-t_0}$$

So your chirp signal is missing a 0.5 correction for the frequency variation.

Use the scipy.signal.chirp:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp

def stft(x, L, overlap, window):
    opoints = int(L*overlap)
    X = np.array([np.abs(np.fft.fft(window*x[i:i+L])[:L//2])
                 for i in range(0,len(x)-L,opoints)])
    return X

tmax = 5
fs = 11025
N = fs*tmax

t = np.linspace(0, tmax, N)
x = chirp(t,0,tmax,2000)
windowsize = 1024

f = np.linspace(0, fs, windowsize)
window = np.hanning(windowsize)
stft = stft(x,windowsize,0.5,window)

plt.imshow(stft.T, cmap = plt.get_cmap('gist_yarg'), origin='lower', aspect='auto',
           interpolation='nearest', extent=(0,t[-1], 0,fs/2))
plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.show()

Sources:

1: https://en.wikipedia.org/wiki/Chirp

2: Basic question about chirp signal


Using only positive frequencies in fourier domain, How will it affect the ifft? I am going to do some kind of transformation and transform a data to another domain, and again back to the first domain. For this, I take a fourier transform of the data and separate the positive part and I use the positive frequencies only.

I know that for example in the frequency domain of a wavelet we have positive and negative frequencies. what will be the effect of using only positive frequencies?

You can do things that way, but unless you're manipulating things properly, you won't get the right answer.  A procedeure for doing a N=32 point example is:

1) generate 32 point time domain data,
2) forward transform it,
3) divide results by 32,
4) multiply positive frequencies by 2 (frequencies f= 1 to 15),
5) zero negative frequencies (frequencies f=17 to 31),
6) inverse transform it,
7) ignore the imaginary outputs.

The real results of the above procedure will match the original input waveform.

Note that in the above, you don't multiply the f=0 or f=N/2 outputs of the FFT by 2.  For this case (N=32 - an even number), they are unique - they don't have negative frequency counterparts.  Now if you're using an appropriate FFT routine that allows for odd N, then you have to account for the fact that you'll have an f=0 point, but no f=N/2 point.

Unless you are doing very large FFT's, you're probably not going to save much in the way of processing time by using just the positive frequencies.


High quality reverb algorithm I did a bit of searching on this site but surprisingly, I didn't find much relevant information, and my knowledge of DSP is very very limited.

My goal is quite simple : I want to program an algorithmic reverb in C++ which sounds really good. Or more exactly, the best option would be to let the end user choose the trade-off between quality and cpu-usage.

From what I found out so far is that, to create a reverb, you must feed the dry signal into a early reflections algorithm, then to a late reflections algo. Is this correct ?

Now, I found an extensive article on the late reflections part, using Feedback Delay Network (REDUCING ARTIFICIAL REVERBERATION REQUIREMENTS
USING TIME-VARIANT FEEDBACK DELAY NETWORKS). From what I've read, FDN is a high quality, not too expansive (cpu wise) way of simulating the late reflexions. Furthermore, I guess you can control the quality/cpu charge trade-off by changing the number of delay lines.

However I have absolutely no idea of how to program the early reflections algo (remember ? I'm  really ignorant in the DSP domain).

It would sound logical to me to use a kind of multi delay, which has the advantage to be easy to program and computationally inexpensive. But that sounds too simple to be true.

Furthermore, my instinct tells me that one or several filters must be included somewhere in the signal path.

Could some one please clarify this topic a bit for me ?

Two notes :

I'm not after a convolution reverb at all. I don't really care about the realism of the reverb, but instead I want a good sounding, tweakable, not to cpu-hungry reverb.
Also, the coding part isn't what worries me, I'd ask on stackoverflow otherwise. It's really the DSP part, and just that part, which I'm after :)

You need to generate early reflections with a few taps of delays (= convolution with the sum of a handful of diracs) ; and the "tail" with what is usually implemented with a network of all-pass (AP) and comb filters.

The first part is trivial to implement but difficult to get to sound right. It might help to look at the positioning of peaks at the head of pre-recorded impulse responses to get an idea of what kind of responses are "natural". It's hard to make it easily parametrizable by a few settings, though you could get away with it by offering a bunch of presets for this part to get an "hybrid" reverb à la Virsyn Reflect.

The tail is where algorithmic reverbs shine (no pun intended), since it is easy to parametrize. The grandfather of algorithmic reverbs is Schroeder's algorithm. Note that it only generates a "tail", not early reflections - you could augment it with a few taps of delays to give it a head. It does not sound very good (very "grainy") but it's a decent start - with that you have something to mess with and it helps understanding the impact of each parameter. Many praised algorithmic reverbs, especially those from the 80s (Lexicon, Eventide, Publison) are just finely tuned topologies of all-pass filters and comb filters. I suspect this involved a lot of trial and error by trained ears who got very good at figuring out how changes in parameters or topology would sound, rather than any scientific approach. Here is another interesting read - showing Keith Barr's favorite reverb topology. The original midiverb might not be your definition of "high quality algorithmic reverb", but it sounded sweet, and it puts thing in perspective to know it was made with a "dsp" only capable of MACs with coefficients of 0.5...

Valhalla room is my favorite software algorithmic reverb - take some time to read its developer's blog to find inspiration about the kind of things that help getting it right.


How can you tell whether there's a signal (when your signal looks a lot like noise)? This is my snoring detector again.
I've gotten pretty good at detecting a signal when there's anything there -- can track from a wall-peeling snore down to breathing you can't even hear in the recording.  The problem is, I can't tell when the signal has dropped below detectable level and the app is just "hearing things".  And, unfortunately, snoring/breathing is often irregular enough that a simple autocorrelation or similar interval timing scheme is unlikely to help much.  (And it's actually likely that in some cases the noise is more regular than the breathing.)
So, are there any tricks I'm missing for figuring out when there is no signal?  It kind of seems that I'm up against a hard place here, given the "signal" is so noise-like to begin with.
(And maybe this is related to another problem I'm having:  Strangely, I can't accurately (or even approximately) measure the signal level even when fairly loud.  Since I need to use rolling averages and ratios to detect the signal anyway, the level information kind of gets lost.  I'm looking for some tricks to reconstitute it.)
Basic technique
(For Yoda)
The audio signal is sampled (generally at 8000Hz, for various reasons), then FFTed in 1024 blocks.  (In my experiments Hamming filters and overlapping blocks seem to have little effect, though those may be revisited later.)
The FFT is divided into "bands" (currently 5, slightly skewed in size to place more detail on the low end) and the "spectral difference" and level of each band is summed.  Long-term averages of the peak-limited values are used as "thresholds", and further bias adjustments are used to maintain a roughly 20% "over threshold" rate.
Each "over threshold" value is given a weight of 1 (under threshold is given a weight of 0), but then that weight is adjusted by the apparent "variability" (at roughly 2Hz) in the band, to give more weight to bands that carry more apparent signal.
The weights of the bands are summed and then the summed weights of subsequent blocks are summed over about a second to produce a running "score".  This is again compared to a running average threshold (plus several heuristics) to detect snore onset/offset.
Update
It suddenly occurred to me that if my algorithm effectively maintains a constant-level signal (per my signal level problem), the way to effectively gauge SNR is by measuring the noise when there's no signal.
Conveniently, snores are intermittent, with lots of "dead air" in-between.  And I'm already detecting the snore envelopes.  So anything outside of the envelope (between the end of one snore and the start of the next) is presumably noise!  This I can (with some modest degree of accuracy/repeatability) measure.  (It took three tries to come up with a halfway decent algorithm, of course -- reality never matches the theory.)
So I don't have the full answer yet, but I've made progress.
(While the above technique gives me a fairly good proxy for SNR, I'm still having trouble estimating actual signal level.  My "relative level" indications can be off the scale for a barely audible breath and so-so for a window rattler.  I need some sort of proxy for absolute level.)

So far you seemed to identify snoring through the presence of a strong periodic pattern in the audio. If you tell me that there might be other sources with the same property, it's time to move on and focus on a property of the signal more specific of snoring ; and I would suggest to look more precisely into the timbre of the recorded sound itself. Given that a human listener can tall apart snoring from breathing, you might resort to traditional audio signal classification approach. Compute features like MFCCs and Fluctuation patterns on your audio data, and train a bayesian classifier (or more sophisticated) to discriminate snoring from breathing.

It might be time also to review the literature on the topic... For example I found this: http://web.itu.edu.tr/kamasak/pubs/pdf/pm7_8_007.pdf


Acquisition of vibration signal and processing help so I've been working on a side project for some time now.

I have rotating shaft that's starting from 15k and goes up to 150k RPM(250-2500Hz), I have the vibration signal of the shaft from transducer and also I have pulse when the shaft completes each rotation (so I can calculate the RPM).

I've been searching the internet for quite some time now reading a lot of theory about dsp and order tracking in general but I simply cannot make it work.

So my questions is, what approach would you use to extract useful information and plot it on a graph using Vibration vs RPM?

(I've been reading something about sampling the vibration signal according to the rotation pulse but I'm just confused is once per revolution enough for the vibration signal, also whats after I have vibration signal and rpm signal)

I will do all the research myself, but i need some help starting out. I would be really happy if someone can write down some general steps I can follow.

I suppose that you are interested in measuring vibration for monitoring purposes, whether for the shaft's bearings or what the shaft itself is doing as it gets loaded (e.g. sagging and causing vibration).

For this purpose, the vibration signal can be assumed stationary. This means that given enough time, no new components will be added to the vibration signal. The vibration is "the same" in every rotation.

With these assumptions in place, what you can do is to integrate an aspect of your vibration signal over short periods of time and cross referencing it with your RPM.

For example, let's say that the vibration signal that you get from your sensor is represented by some $x(n)$. If you want an estimate of how strong the vibrations get, you can simply integrate the $|x(n)|$ over time (that is, the absolute value of $x$) and see if the amplitude of the vibration seems to increase with respect to RPM. The integration itself can be done with a simple moving average filter. If your Sampling Frequency is 7kHz ($F_s > \frac{f_h}{2}$) and your highest frequency $f_h$ here seems to be 2.5kHz), then this means that 1 revolution at 250Hz takes about 140 samples or approximately 20ms. So, if you were to filter your $|x(n)|$ with a simple moving average filter with length 140 samples, you would generate another signal that would give you average amplitude of vibrations per revolution at your low end RPM. As the RPM increases of course, the smoothing will be greater.

But, you are also very lucky because in addition to your vibration signal you also have the RPM sensor. From these two and the point of view of monitoring the bearings you can understand two things: Whether the bearing's shell has developed a chip, or whether one of the bearing elements is deformed. Here is why: If the bearing shell has a chip then there will be a peak in the vibration signal at the same phase in every revolution. If one of the bearing elements is defored (whether cylinder, ball, other) then, as the bearing rotates, the elements meet the shell at a different point at each revolution. Therefore, your get a "rolling" pulse in the signal. On the first revolution it occurs 2ms after the rotation starts, on the second revolution it occurs 4ms after the rotation starts, on the third revolution it occurs 6ms after the rotation starts and so on.

To see these things, you can plot your vibration signal in parallel (or, in phase) with your rotation signal and look for "regular bumps" or "rolling bumps".

The same applies for when you want to measure if your axle is sagging beyond an RPM (because it was asymetrically lathed or it was loaded suddenly or whatever). You just monitor the vibration signal for bumps versus the RPM.

Towards that direction, you might find wavelet analysis more useful.

Hope this helps.


library for computing differential signatures from 2D planar curves For one of the lectures on 2D planar differential geometry in the series "Image and Video Processing: From Mars to Hollywood With a Stop at The Hospital", the lecturer discusses differential signatures based on parameterization of a 2D planar curve (starts at approx minute 17 in the link attached). I'm not sure if there is some other terminology used but I'm having a hard time finding a software library that can extract such differential signatures once a curve has been segmented/identified in an image. Is there some other name that such a process is called or is there some other method that can extract similar features? I tried searching the OpenCV library but I wasn't certain.

Recently I have worked on a projective invariant line signature:

  Projective splines and estimators for planar curves
  Thomas Lewiner and Marcos Craizer International Journal of Computer Vision, 2000
  http://thomas.lewiner.org/pdfs/projective_splines_jmiv.pdf

That paper proposes a projective length and a projective curvature estimators for plane curves, when the curves are represented by points together with their tangent directions. The goal of the paper is an invariant spiral construction in the canonical projective frame.

If you like to implement this paper, the paper is pretty self contained. And I posted the details which are kind of hidden in my blog along with some code.

If you don't need projective invariance then I would suggest Fourier Descriptors or Elliptical Fourier Descriptors. There are many MATLAB implementations available for that.

And this work might be quite nice to take a look, if you are concerned about numerics and invariance.


Bit reverse order Technique in FFT Can anyone tell me how does bit reversal order technique works for FFT when you are trying to break down the FFT to small sizes. Like I want to only for the when the index is odd because that involve some reverse carry propagation, but when the index is even you just add N/2 to it to get your reverse index r. An example of N=8 would be perfect if someone can explain me.

If I understand your question correctly you want this (index | binary | bit rev. | bit rev. index):

$$
0\quad        000\quad      000\quad 0\\
1\quad        001\quad      100\quad 4\\
2\quad        010\quad      010\quad 2\\
3\quad        011\quad      110\quad 6\\
4\quad        100\quad      001\quad 1\\
5\quad        101\quad      101\quad 5\\
6\quad        110\quad      011\quad 3\\
7\quad        111\quad      111\quad 7\\
$$


How can I transform a time signal into a third octave spectrum? I try to make a third-octave spectrum from a time signal and I have some troubles with the algorithm. I'm not an expert in signal processing so I don't know what is wrong.

Here what my algorithm looks like :

Time signal -> Z transform of the signal -> First octave band filtering -> Quadratic sum of the result's values => First octave band value.
Under factor 2 sampling ->Z transform of the signal -> Second octave band filtering -> Quadratic sum of the result's values =>Second octave band.

Now I have this, I try to test it on Matlab with a time signal about 2 millions points.

I have this code to calculate all the octave's bands :

    fc = fMin*2.^((0:bmax)*bw);             %Center frequencies.
    fl = fc*2^(-bw/2);                      %Lower cutoffs.
    fu = fc*2^(+bw/2);                      %Upper cutoffs.

And this one to calculate the Z transform of the time signal :

    x = [1 2 3];
    syms z;
    d = z.^(-1*(0:numel(x)-1));
    zTransform = sum(x.*d)

My main problem is, I'm not completely sure that my algorithm is right and if it is, how can I filter the z transform equation that the code above gives me with each third-octave band filters?

Thanks in advance, it's been a while I'm on this problem.

maybe either the bandpass filters or the peaking EQ filters in the audio EQ cookbook can help.  bandwidth is defined in octaves in that.  the output, $y$, of each 1/3 octave BPF filter should be squared and low-pass filtered to get the mean square, $\overline{y^2}$.  and then

$$dB = 10 \log_{10}\Big(\,\overline{y^2}\,\Big)$$

is applied to that to get dB.

are you planning on doing this to the whole audible spectrum?  31 bands?


Image downsampling with Haar DWT? How can the resolution of an image in matlab be reduced by using 2D haar wavelet by downsampling-decimation? Is there any predefined code or tool within wavemenu? thanks in advance

yes ok thanks for the reply. Specifically, what I intend to do is an algorithm of digital image processing for use in retinal implants. Clarifying a bit my problem I know that there is a mathematical operation applicable to digital images that is called Wavelet Transform and that is in some way to compress the image can also lower its resolution through downsampling. My question is how can I from Matlab (with wavemenu GUI for example) compress an image using the Wavelet Haar transform and lower its resolution for example 640x480 to 32x32. In the following article is explained in a general way in the last step of table 1
https://www.researchgate.net/publication/323193463_Comparative_analysis_of_image_processing_algorithms_for_visual_prosthesis