Why is the gain of my IIR filter positive? Sorry, two questions in one day!

I'm struggling to understand what I'm doing wrong in this very simple filter design. I want to design a simple, single pole low pass filter and implement it as an IIR in an embedded system.

I was looking at the DSP book online at recursive filters, and using their single-pole example to get started: http://www.dspguide.com/ch19/1.htm

According to the next page the coefficients I should need for the IIR single pole low pass filter are:
$a_0=1-x$ and $b_1=x$

In order to get $x$ they give the equation (19-5) of $x=\epsilon^{-2\pi f_c}$

Given an fc of 0.125, I calculate x = 0.455938. Is this correct?

An fc of 0.125 at a sampling rate of 8kHz is 1kHz.

I don't see why when I try to get the frequency response of this filter I end up with entirely positive gain, as shown below:

http://www.valvers.com/wp-content/uploads/2013/01/butterworth-lowpass-positive-gain.jpg

Here is the code I've been using:

%Start from nothing!
clear;

% Set the sampling frequency used by our digital filtering system:
fs=8000;

% Set the coefficients up (As already worked out!) for a single-pole low pass
% filter. This should give us -3dB @ 1kHz with -6dB/Octave roll off
a = [ 0.544062 ];
b = [ 1, 0.45594 ];

% Determine the frequency response of the filter design above. Get the output
% in frequency rather than rad/s. Use 64 plot points.
[H,f] = freqz(b, a, 64, fs);

% Show our cutoff frequency
cutoff = -3 * ones(64);

% Plot the result so that we can see if it is correct:
figure(1);
plot(f, 20*log10(abs(H)), f, cutoff);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');

I'm again, just stumped as to what's going on, and don't know how to solve the issue.

You have your a and b vectors reversed. For a first-order IIR filter, the a vector should have two elements (and you almost always want the first element to be unity). So, you're plotting the frequency response of the filter with difference equation $0.544062y[n]=x[n]+0.45594x[n−1]$, which is not what you want.

Also, as a note on MATLAB syntax, you'll also need to negate the second element of a (multiply it by -1) in order to get what you want. Read the documentation on functions that take IIR filter coefficients to learn why.


data for panorama stitching testing I want some data (images that have only shift) to test on xuvtools  or Microsoft ICE and my own software.
Is there any available?

From a simple Google search...

http://sourceforge.net/adobe/adobedatasets/panoramas/home/Home/
http://www.ti.uni-bielefeld.de/html/research/avardy/index.html
http://sun360.mit.edu/
http://www.iiia.csic.es/~aramisa/datasets/iiiapanos.html
http://www.visualsize.com/mosaic3d/index.php


Signals - Removing Blocks I'm after some advice on a problem that I am facing regarding signal processing.

I am processing a voice clip (.wav) and this voice clip contains noise, and I need to remove the noise and just concentrate on the Phones so I can analise them.

What I have thought about, and, what I have researched is, if I use zero-crossing as well as the total energy I can remove the blocks that do not contain sufficient energy or frequency.

The problem at the minute is that I have a 1D vector, acting as a 2D vector and let's assume the following data is used:

v1 = {128, 128, 128, 127}

Now my current algorithm squares each of these numbers and then sums the entire vector up. This, would therefore not be able to give me an accurate identification of noisy signals. So should I therefore split the signal into blocks (In that I get a 2D vector)? e.g.

B1 = {128, 128};
B2 = {128, 127};

e1 = {sum( (128^2) (128^2) )};
e2 = {sum( (128*2) (127^2) )};

Therefore, this would mean I am given 2 (In this example) energies which I can then check to see if they match the threshold.

Last question, do you think the methods of zero-crossing AND energy are good an accurate way of determining noise and removing it?

Hope someone can help me, thanks :)!

An energy or zero-crossing criterion might help in detecting segments containing no voice, but it won't help in removing the noise superimposed with speech. Since you are mentioning removing noise for improving speech analysis, this doesn't look like the right technique to use.

If you compute your energy using small windows, down to summing a few adjacent coefficient, you will see fluctuation in energy which are just due to the periodic variation of the waveform itself! As a rule of thumb, the window size on which you compute an energy feature must be at least one magnitude bigger than the period of the signal to analyze. Assuming a f0 of 200 Hz, this gives you a minimal integration window of 50ms, or 400 samples at 8kHz.


Apply decay envelope to waveform without multiplication I'm trying to apply a decay envelope to a sinusoidal waveform using Verilog. Hardware constraints prevent use of multiplication to simply multiply by the envelope. The sine values and envelope values, for which I'm simply using $1-e^{-x}$, are stored in ROMs to speed everything up.

Right now I'm simply calculating the output sample by subtracting the exponential value from the sine value if the sine value is larger. If the exponential value is larger than the sine value, the output sample is zero. This works pretty well, but the current exponential value to subtract needs to be modulated by the current location on the sine wave so that the full value is subtracted when the sine is maximum and zero is subtracted when the sine is zero.

I'd appreciate any advice that anyone can give and I can provide more details if required.

There are other methods if you can't support direct hardware multiplication, but you're going to have to sacrifice something. For instance, if you have multiple clock cycles to spare, you could use an iterative approach. C-like pseudocode:

accum = 0
for (i = 0; i < LENGTH(envelope); ++i)
{
    if (((envelope >> i) & 0x1) == 0) accum += signal >> (i + 1)
}
// "accum" now holds the multiplication result

Translating the above functionality to Verilog is an exercise left to the interested reader (it's been a while since I've written in that language). This scheme essentially performs long base-2 multiplication using additions and shifts. For each bit that is set in the envelope scaling factor, an appropriately-right-shifted version of the signal sample is added to the accumulator. The result of the multiplication is held in the accumulator after the loop is complete.

As I said, this isn't the most efficient approach in the world; it's probable that a good synthesis tool would come up with something better even on a device without hardware multipliers. This method trades resource usage for latency in the form of extra clock cycles required to perform the iterations of the unrolled multiply. I'm left to wonder what drives the no-multiplication constraint in your application. Could you add more details on that?

Note: You may need to tweak the amount of extra shift that you apply to the signal when adding to the accumulator based on the numeric format of accum and signal. If they are both of the same fixed-point format (like Q15), then the above approach should work properly.


Envelope follower not smooth enough With the following code, I get the "envelope" of a sound signal.
However in the "envelope" there are small oscillations (high frequency) : it is not smooth enough...

I imagine that it's possible to apply a low-pass to the "envelope signal" in order to  remove the high frequencies oscillations, but I also imagine it's probably better to do it directly by modifying this code :

double attack_coef = Math.Exp(Math.Log(0.01)/( 10 * 44100 * 0.001));   // 10 ms
double release_coef = Math.Exp(Math.Log(0.01)/( 100 * 44100 * 0.001));  // 100 ms
double envelope = 0.0;

for (int n = 0; n < file.Length; ++n)
{
 double tmp = Math.Abs(x[n]);
 if (tmp > envelope)
  envelope = attack_coef * (envelope - tmp) + tmp;
 else
  envelope = release_coef * (envelope - tmp) + tmp;

 env[n] = envelope;
}

Do you have an idea for that ?

You are already just low-pass filtering the magnitude of the input signal. Right now, your envelope detection "circuit" is just a single-pole filter. Let $f[n]$ be the input to your system and, for simplicity sake, assume that attack_coef $=$ release_coef $ = a_1$. Then we can express the output $y$ of your system as follows:

$$y[n] = a_1 \left(y[n-1] - |f[n]|\ \right) + |f[n]|$$

Which we can rearrange as a transfer function $h$ as follows:
$$y = a_1 \left(y\,z^{-1} - |f|\right) + f$$
$$h = \frac{y}{|f|} = \frac{1 - a_1}{1 - a_1 z^{-1}}$$

You can use this information to either adjust your coefficients to get the desired cutoff frequency, or you could add another pole if you want a steeper roll-off.


State of the art motion blur recovery Could the experts point me towards the best motion de blur algorithm? My images are blurred where I know the motion information with good precision(direction and angle)

There are so many deconvolution algorithms however I was not able to figure out which one would be best.

It's really hard to say which one is the best unless we define more precisely what we mean by "best." Some algorithms measure performance using MSE, others use more fancy metrics like entropy.

Since you know the angle and direction of blur, it seems you have a good idea about the blurring function (point spread function). I would start off with a minimum mean squared error type filter (Wiener filter).


Real-time block processing setup on ARM Newbie question here.

I'm working on a system that is running a 32khz sampling rate, audio, 16 bit DAC.
Previously, our system had a very hard latency requirement of 1 sample. So no problem, we just read 1 sample in, did some processing within the sampling rate, and then read 1 sample out.

Now that this requirement has been relaxed, I think that we can use block processing to do things such as spectrum computation and processing and use more than 1 sample to do the adding and multiplication, etc.
Does this buy us blockSize * (1/fs) compute time?

My problem is, being a noob, I don't know how to formulate this new requirement in our chip's main() function.

So let's pretend we have some function SPI_Read(&mic) which is reading a sample, and SPI_Write(&output) which sends to a speaker. How would I transform this into a block based system?

The current system is setup like:

main() {
    while() { // check for interrupt
        // while loop executing on system timer set to 1 / fs
        SPI_Read(&in); // read from adc
        // filters
        SPI_Write(&out); // output to dac

    }
}

This is running on a cortex-M4, no codec.

Thank you for your help! (:

First, use a timer and an ISR to get accurate timing (don't forget to configure the NVIC so that this timer interrupt takes over any other ISR that would be running). Only this will ensure a consistant sample rate. Little variations in timing would create noticeable degradations of audio quality. In particular, in your current example, unless the "filters" code has no branches, the time at which the SPI_Write occurs will be jittery depending on the path taken in the code.

For block processing, you need a "ping-pong" buffering scheme.

In your per-sample ISR:

SPI_Write(output[active_buffer][sample_index]);
SPI_Read(input[active_buffer][sample_index]);
++sample_index;
if (sample_index >= BUFFER_SIZE) {
  active_buffer = (active_buffer + 1) % NUM_BUFFERS;
  sample_index = 0;
}

In your main code:

while (render_buffer != active_buffer) {
  // Process samples from input[render_buffer][0] to input[render_buffer][BUFFER_SIZE - 1]
  // And write them to output[render_buffer][0] to output[render_buffer][BUFFER_SIZE - 1]
  render_buffer = (render_buffer + 1) % NUM_BUFFERS;
}

render_buffer and active_buffer are declared as volatile. Start with NUM_BUFFERS equal to 2. Higher values will introduce more latency, but might be necessary if the processing code has exceptionally long code paths.

The benefit of processing samples block by block is that some variables required by your DSP algorithm (such as filter state variables) can be kept in registers within the processing loop, saving load/stores. There are also situations in which you can benefit from batch load/stores. Finally, it makes it easier to schedule computations at "audio rate" (once per audio sample) and "control rate" (once per block).


Derive minimum phase from magnitude With the desired magnitude of a transfer function in the frequency domain in C++ as described below what is the correct corresponding minimum phase?  In general how does one derive the correct minimum phase in the frequency domain if the desired magnitude is known, specifically such that the a complex IFFT of it would return an entirely real part.

// PnCutoff.cpp : Pink noise above a cutoff.
#include <tchar.h>
#include <math.h>
#include <conio.h>

class PnAtFilter
{
    double _fcHz;
    public:
    PnAtFilter(double fcHz) { _fcHz = fcHz; }
    double Mag(double fHz)
    {
        if(fHz < _fcHz) return 1.0;
        else return  1.0 / sqrt(fHz/_fcHz);
    }
    double Phase(double fHz)
    {
        return 0.0; // what is the the correct minimum phase in r?
    }
};

int _tmain(int argc, _TCHAR* argv[])
{
    PnAtFilter* pnf = new PnAtFilter(1000);
    for(int fHz = 10; fHz <= 100000; fHz *= 10)
    {
        double mag = pnf->Mag(fHz);
        double phase = pnf->Phase(fHz);
        _cprintf_s("%d Mag=%f, Phase=%f\r\n", fHz, mag, phase);
    }
    _getch();
    return 0;
}

/* output:
10 Mag=1.000000, Phase=0.000000
100 Mag=1.000000, Phase=0.000000
1000 Mag=1.000000, Phase=0.000000
10000 Mag=0.316228, Phase=0.000000
100000 Mag=0.100000, Phase=0.000000
*/

Check out Julius O. Smith III's write up.

There is a Hilbert transform relationship between the magnitude response, $G(\omega)$, and the phase response, $\theta(\omega),$ of the associated minimum phase filter.

If

$$
H(\omega) = G(\omega) \exp(\jmath \theta(\omega))
$$
then
$$
\ln( H(\omega) ) = \ln(G(\omega)) + j\theta(\omega)
$$
and
$$
\theta(\omega) = - {\cal H}\left[ \ln(G(\omega)) \right]
$$
where ${\cal H}$ is the Hilbert transform.

Also check out Wikipedia's write up.


Matlab: Calculation range I have this Matlab/Octave code:
% needed for Octave ----------------------
pkg load signal
% ------------------------------------------------

% Square -----------------------------------------
figure
N=50;                               % harmonics
fs = 1000;
r_start=0;
r_end=10;
r_step=1/fs;
r = r_start:r_step:r_end;           % range
w_sqr = square(r)/2;                % Square wave [-0.5:0.5]
%w_sqr = square(r);                  % [-1:1]

plot(r,w_sqr);
axis([r_start r_end -1.2 1.2]);
grid on;
hold on;

% Fourier -----------------------------------------
i=1;
sum=0;
for t=r
  for n=1:N
    sum = sum + (2*sin(n*t)+sin(pi*n-n*t)-sin(n*t+pi*n))/(2*pi*n);
  end
  F(i)=sum;
  i=i+1;
  sum = 0;
end;
F=F';
plot(r,F);
axis([r_start r_end -1.2 1.2]);

which results this:

By changing the value of variable w_sqr to (square(r)+1)/2 and calculation of F to F(i)=(1/2)+sum I get the result in range [0:1].
What changes are needed to do in Stage 2 to get the result in range [-1:1]?


How is cross-correlation related with orthogonality? In linear prediction we can say that in case of optimum linear predictor the error with be orthogonal to data. And when we derive minimum mean square error for $\underline{y} = \mathbf{a}\underline{x} + \underline{b}$ case we find following relation

  $\xi = \sigma^2_y (1 - \rho_{xy} )^2$

Where $\rho_{xy} $ is correlation coefficient and if we talk about cross correlation $r_{xy}(k,l)$ of zero mean and uncorrelated random variables x and y, then how does this conclude that if $r_{xy}(k,l)$ will be then the signals will be orthogonal?

I suppose you mean the cross-correlation at lag zero. Well take an Hilbert space $H$ (i.e. a metric space in which you can define a scalar product $\langle\cdot ,\cdot\rangle$). Then $x,y\in H$ are orthogonal if $\langle x,y\rangle=0$, by definition.
If your Hilbert Space is $L_2(\mathbb{R})$ (the space of real square integrable functions) then the scalar product of $f,g\in L_2(\mathbb{R})$ is defined as
$$
\langle f, g\rangle = \int_{-\infty}^\infty \overline{f(t)}g(t) dt
$$
If the signals are assumed to be ergodic, which is actually a costumary assumption, then the cross-correlation is simply a sliding scalar product, you have
$$
\text{corr}(f,g)(\tau) = \int_{-\infty}^\infty\overline{f(t)}g(t+\tau) dt
$$
hence
$$
\text{corr}(f,g)(0) = \langle f ,g \rangle
$$
and uncorrelation means orthogonality.
The same applies for discrete time signals in $l_2$


How do I compute the impulse response of a robot? I am trying to model a robot (my question in robotics.SE), so that I may train it using some reinforcement learning techniques. I was wondering if I can compute an impulse response of a four wheel robotic chassis empirically? Or something similar that might help me model the robot. Does a thing like 'the impulse response of a robot' even make sense?

Since a single strobe of a PWM signal shall hardly cause any movement, but a steady stream of PWM signals will, can I consider the impulse response as an infinite impulse response that superposition to cause perceivable motion?

I am certain that the robot chassis is a non-linear system. Generally speaking, will impulse response make any sense on a non-linear system?

Impulse response makes sense for linear systems. If your system is nonlinear forget it. If you assume that your system is linear time invariant, you may want to measure the impulse response by exerting a short duration high energy pulse to your robot. It must have sufficiently high energy so that most of the exerted energy is transferred to kinematic energy of the robot to enable motion.

If this is not possible you can apply longer duration pulses to your robot and then measure its response (motion). You must deconvolve then the response and the input to get the impulse response of your robot. You may try this for a number of different shape pulses and if you get impulse responses similar to each other every time, you can be fairly sure that your system is linear.

Your comment on having an infinite impulse response is not in line with correct systems modelling, so forget about it.


DFT coefficients meaning? What "are" they? What's a sensible way to interpret the coefficients (and what isn't)? To pose specifics:

DFT coefficients describe the frequencies present in a signal
They describe the sinusoidal frequencies of the source

Are either of the above true?

At it's most fundamental, the DFT is about fitting a set of basis functions to a given set of sampled data.  The basis functions are all sinusoidal functions, expressed as the complex exponential with a purely imaginary exponent.  Using the most common scaling convention each basis function, without its scaling coefficient, is:
$$ g_k[n] \triangleq \tfrac1N e^{j \omega_k n} $$
For the DFT:
$$ \omega_k = \frac{2 \pi k}{N} $$
and the basis functions add up as
$$\begin{align}
 x[n] &= \sum\limits_{k=0}^{N-1} X[k] \ g_k[n]\\
\\
  &= \sum\limits_{k=0}^{N-1} X[k] \ \tfrac1N e^{j \omega_k n} \\
\\
  &=  \tfrac1N \sum\limits_{k=0}^{N-1} X[k] \ e^{j \frac{2 \pi k}{N} n} \\
\end{align}$$
It's real easy to solve for the coefficients:
$$\begin{align}
 X[k] &= \sum\limits_{n=0}^{N-1} x[n] \ e^{ -j \frac{2 \pi k}{N} n} \\
\end{align}$$


Frequency of a noisy signal with fixed width pulses My electricity meter has a front panel LED which pulses in proportional to the current energy rate. The frequency changes but the pulse width is set at 250ms.  Due to environmental conditions there's plenty of noise. Currently I'm sampling this at 10kHz, the pulses may be a few seconds apart.

I would like to filter these pulses and ignore anything else that wasn't 250ms long (ie the noise) and find the frequency of the train. Are there any suitable methods to achieve this?

My first thought was to differentiate the signal but this would be easily confused by large noise spikes.
I also wondered about testing it's coherence with a simulated pulse train, but as the frequency changes this would be potentially infinitely intensive AFAIK.

If I was to build an electronic solution I'd use edge comparators to find the pulses.  I'm wanting to use NumPy/matplotlib and deal with the signal in the frequency domain, if possible (purely as an educational exercise).

Not an easy problem to solve...

You could always cross-correlate the signal with a replica square wave and then apply simple thresholding to find the pulse locations. Then just count the peaks in a given interval to determine frequency.

I think the square wave nature of the signal could make the signal quite difficult to deal with in the frequency domain because of all the harmonics that comprise it.


Signal Processing using Fourier Transform So I'm trying to understand how MRI machines work. I understand all the concepts of it, the parts, what they do, how the machine works, etc. The part I'm having trouble with is the fourier transform business.

In the MRI machine, when the hydrogen atoms release the radio frequencies previously absorbed, the transceiver absorbs it and uses fourier transform to convert those signals from the time domain to the frequency domain.

The part I don't understand is, how does the fourier transform end up producing the image of the body. How do you go from math to a detailed picture of the brain??? Can anyone explain that?

The signals that we measure in MRI are a combination of signals from all over the object being imaged. It so happens that any signal (even if you simply make one up and draw a squiggle) is composed of a series of sine waves, each with an individual frequency and amplitude. The Fourier transform allows us to work out what those frequencies and amplitudes are. (That is to say, it converts the signal from the time domain into the frequency domain.)

Since we encode the signal with magnetic field gradients which make frequency and (rate of change of) phase relate to position, if we can separate out the frequencies we can say where we should plot the amplitudes on the image.

The amplitudes serves as the brightness levels in the image:

$Amplitude = I$

$(Phase, Frequency) = (x,y)$

$I(x,y) = I(x i + y j)$

Phase and Frequency help to find where in the image(which exact part of the  MRI) has the particular amplitude measured by the technique of passing RF signals through the body and collecting the signal information.

The exact process is well explained with animations at:
http://www.imaios.com/en/e-Courses/e-MRI/Signal-spatial-encoding

An RF pulse does not have one frequency only (for this, it would need to be of infinite duration). It covers a certain bandwith, which depends on the shape of the pulse and its duration.

The thickness of the slice can be varied by adjusting the bandwidth of the selective pulse  and the amplitude of the slice selection gradient.

For a fixed amplitude gradient, the wider the bandwidth, the greater the number of protons excited and the thicker the slice
For a fixed bandwidth, the stronger the gradient, the greater the variation of precession frequency in space and the thinner the slice

Moreover, the shape of the RF pulse in time will also determine the bandwidth profile in frequency, and thus the slice profile.


FFT: Extracted spectrum depends on number of points sampled I have a signal $a(t)$, this is how it looks when I plot it:

Zoom in

I am trying to extract the eigenfrequencies of this signal using FFT in MATLAB. I do it like this

Y=fft( a(1:3:end) );

Fs = 1/1.5;       % Sampling frequency
T = 1/Fs;         % Sampling period
L = dms(2);       % Length of signal
t = (0:L-1)*T;    % Time vector

P2 = abs(Y/L);
P1 = P2(1:floor(L/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
plot(f, P1)

Now, it turns out that if I use Y=fft( a ), I only get frequencies close to 0 when I plot the spectrum:

If I use Y=fft( a(1:3:end) ), then I get frequencies at 0 and nonzero frequencies:

Why do these nonzero frequencies not appear when I do the FFT on the full dataset, but only when I sample every 4th point?
Is the nonzero frequency actually contained in my dataset or is it an artifact?

If I use Y=fft( a(1:3:end) ), then I get frequencies at 0 and nonzero frequencies:

Well, yes, but that means you decimate your signal by a factor of 3, which on one hand will shift frequencies to $\frac13$ of their original frequency, but also will lead to aliasing. In other words: unless you know there's no signal above $f=\frac16f_\text{sample}$, you simply can't do that.

EDIT (based on zoomed-in plots):

So, hand-counted from your most-zoomed in plot, it seems there's about 5 peaks in 20 samples, which implies a period of 4 samples.
That is ok by Nyquist's sampling theorem.

Now, if you just look at every third sample, which you do by doing a(1:3:end), then there's only ~1.333 samples per period. Nyquist demands that there's at least 2 samples per period. So, your second attempt aliases the signal from the $\frac15f_\text{sample/2}$ onto a different point in the spectrum, where it actually isn't. You damaged your signal by undersampling; you mustn't do that. It's not very clear to me how the idea to do that happened in the first place, but you should definitely read up on the Sampling Theorem.


Note Recognition Software I'm wondering if anyone knows of any open source solutions to note recognition. I'm planning on building a program for online sheet music, which recognizes your notes, and tries to tell you what you can improve, and how.

However, for piano for example, I don't even know if it's possible to have software that detects various notes at once, with very high speed. If anyone knows anything like this, please let me know!

Although it looks like it's not being actively maintained, you might check out "aubio" to see if it will meet your needs : http://aubio.org/

  aubio is a tool designed for the extraction of annotations from audio signals. Its features include segmenting a sound file before each of its attacks, performing pitch detection, tapping the beat and producing midi streams from live audio.


Tune a piece of music Suppose I have an instrument which is out of tune. For example lets say, as an example, the notes it plays are as follows: C+20 cents, C#+20 cents, D+20 cents D#+20cents etc. In other words, it plays all its notes 20 cents too high.

I then record this instrument, having played polyphonic music on it.

What algorithm can I use to tune the track so that the notes become: C+0 cents, C#+0 cents, D+0 cents, D#+0 cents. I wish to solve this in the frequency domain and I have implemented an FFT and IFFT in this regard.

UPDATE

Please note that the algorithm should work out the 20 cent offset by itself without me having to supply it. I might later record another instrument with a -35 cent offset. All I know is that the notes are relatively in tune, but the instrument itself relative to other instruments might be out of tune by x cents.

I am interested in human voice, but also other instruments like violin, flute, guitar, piano and analog synth.

Essentially what you'll need to do it adjust the pitch of the audio data such that your slightly sharp notes are 20 cents flatter.

By far the easiest way to do this would be to play the audio at a lower sampling rate, which will lower the pitch accordingly.  However, this approach will also slow down the original tempo.

A solution in which the signal retains its speed would be to use a time-scale modification algorithm to speed up the audio while retaining the pitch, followed by playing at a lower sampling rate, which returns the audio to the original tempo, with the desired lower pitch.

Clearly, the factor by which you adjust the time scale and pitch adjustment factor must be carefully chosen to yield the desired results.

In terms of implementation, playing back a reduced sampling frequency is (relatively) trivial.  A frequency domain algorithm which can perform the time-scaling for you is the phase vocoder.  You might check out the following links to help you implement it.

http://homes.esat.kuleuven.be/~bdefraen/Laroche.pdf
http://www.ee.columbia.edu/ln/labrosa/matlab/pvoc/

Good luck!


Sine of frequency 0 contains sines of all frequencies at once in it You know that a sine corresponds to a pulse by J.Fourier transform. The lower is the frequency, the closer is the pulse to the origin. A constant signal is a sine (or cosine, that may be important) of frequency 0. It is a pulse in the origin. This is ok, since we know very well that a pulse transforms to the white spectrum: it is a combination of all sine waves at once.

But, I see a contradiction here. I have just shown that a constant is also a sine of frequency 0 and, also, a combination of all frequencies. A sine of frequency 0 is a sum of all frequencies. How is this possible?

You are confusing time domain and frequency domain. A constant time domain function (a sine with frequency 0, if you like) corresponds to a delta impulse at the origin (i.e., frequency zero!) in the frequency domain. A delta impulse at frequency zero is zero for all other frequencies $\omega\neq 0$. Consequently, a constant in the time domain is not a combination of all frequencies. What you were probably thinking of is a delta impulse in the time domain, which corresponds to a constant in the frequency domain (i.e., it contains all frequencies).

So, summarizing, a constant in one domain corresponds to a delta impulse in the other domain. A constant in the time domain does not contain any frequencies other than 0 (delta at zero in the frequency domain). A constant in the frequency domain (i.e. a combination of all frequencies) corresponds to a delta impulse in the time domain.


How do you calculate/compute/determine the spectral coverage of a frequency spectrum? I'm comparing frequency spectrum of FSK, ASK, QPSK, and MSK. And I'm asked which has the least and greatest spectral coverage. Could someone enlighten me on how to get it?

Spectral coverage, I assume, means the Bandwidth required to transmit the incoming signal. A simple formula(approximate) for the Spectral Efficiency of Memoryless Modulation Schemes(MSK is not one) which gives the ratio between Data Rate and Bandwidth required is:

R/W = (2*logM)/N

where R is the incoming Data Rate and W is the required Bandwidth; log is in Base 2, M is the number of points on the Signal Constellation and N is the number of dimensions required to use the particular signalling.

FSK requires as many as Dimensions as the number of Points you want to transmit.
ASK is one-dimensional.
QPSK is two-dimensional.
MSK is more spectrally efficient than the other 3 and I think the formula for it's Power Spectral Density and its 3dB value can be found in all Digital Communication Texts.


Kalman Filter Sensor Processing Kalman filter achieves convergence of state vector by using sensor observations.

Assuming a sensor such like velocity sensor, giving two axis velocity information in X-axis as well as Y-axis(see edit for one more example to capture essence of my question)

I can always design kalman filter in such a way that X-axis and Y-axis velocity are taken as two different scalar sensor readings and also a kalman filter where I take them as a vector sensor reading.

From both theoretical and practical aspects Is there any thing I am losing by taking X-axis and Y-axis as independent readings?

Or taking them together helps in converging faster?
What is a better approach? And Why?

Edit: Assume one is in a airplane and a camera captures an image and tells latitude and longitude of nadir. With this data I would like to correct my State vector. Now if I use latitude as a different reading and longitude as a different reading than using both of the data. Will I lose anything? As pointed by an answer below if the signal model has an interaction among the components of reading I lose that interaction. Now do latitude and longitude have an interaction? How do I capture this interaction mathematically with H matrix?

As with any question like this, the answer is: It depends.

What does it depend on? Your signal model.

If your signal model generates $X$ and $Y$ axis velocities independently from each other so that there is no transfer function between the two channels, then two 1D Kalman filters will work the same as one 2D Kalman filter.

If your signal model allows interaction between $X$ and $Y$ axis values, then you lose that captured interaction between the axes, which can degrade the performance of the two 1D Kalman filters over a single more accurate 2D Kalman filter.

To give any more detailed answer will require you to state your signal model.

The edit really doesn't give me enough information, but let me try to update my answer.

I'm assuming, from what you say, that the latitude and longitude are estimated from a nadir image.  If that's true, then unless the image pixels are precisely aligned with the lines of latitude and longitude, you will get cross-coupling between your latitude and longitude estimates.

That cross-coupling means you'll be better off using a 2D Kalman filter.


Do energy of signal and power spectrum (in db/hz) indicate towards same signal property? If I calculate energy of signal by squaring the magnitude of signal sample and power spectrum using spectrogram (db/hz) what kind of property of signal do they indicate and are they point towards same characterstic. That is higher energy signal will always have higher peak spectral power and vice versa. Also what is there practical significance.

On a straight forward approach, energy of a signal given by $E = \int \limits_{-\infty}^{+\infty}|x^2(t)|dt$ is a constant number whereas PSD can be obtained from FFT which is a function of frequency. PSD is the square of the absolute value of FFT.

In either case the value of energy or PSD gives an insight of the amplitude/strength of a signal.


90 Degree Phase Shift / Rotation Algorithm for SDR I'm working on a library to preform frequency demodulation from an SDR signal. I'm using https://github.com/sysrun/rtl-sdr/blob/master/src/rtl_fm.c as a crib for this. The algorithm starts by taking the samples (I & Q data in as 8-bit unsigned values I believe) from the radio and putting them through a 90-degree rotation. From this answer to a related question, I believe I understand why this is needed, but I don't understand how the algorithm achieves this. The algorithm is as follows, where buf holds the (I & Q) samples directly from the radio:

void rotate_90(unsigned char *buf, uint32_t len)
/* 90 rotation is 1+0j, 0+1j, -1+0j, 0-1j
   or [0, 1, -3, 2, -4, -5, 7, -6] */
{
    uint32_t i;
    unsigned char tmp;
    for (i=0; i<len; i+=8) {
        /* uint8_t negation = 255 - x */
        tmp = 255 - buf[i+3];
        buf[i+3] = buf[i+2];
        buf[i+2] = tmp;

        buf[i+4] = 255 - buf[i+4];
        buf[i+5] = 255 - buf[i+5];

        tmp = 255 - buf[i+6];
        buf[i+6] = buf[i+7];
        buf[i+7] = tmp;
    }
}

The comment is somewhat confusing, as I don't understand how the complex values match the array. I assume the values in buf correspond to the parts of the complex numbers like such: buf[i] + buf[i+1]j. If that's the case, then why is this a 90-degree rotation? Why does it require 4 complex values?

(buf[i+0] + buf[i+1]j) -> ( buf[i+0] + buf[i+1]j)
(buf[i+2] + buf[i+3]j) -> (-buf[i+3] + buf[i+2]j)
(buf[i+4] + buf[i+5]j) -> (-buf[i+4] - buf[i+5]j)
(buf[i+6] + buf[i+7]j) -> ( buf[i+7] - buf[i+6]j)

Any insight into this algorithm would be greatly appreciated... thanks!

I believe this algorithm isn’t a rotation (phase shift) by 90degrees but instead is a frequency shift of fsamp/4.
In this algorithm each sample is rotated 90 degrees with respect to the last one causing a frequency shift:

Samp1: no rotation
Samp2: 90 degree rotation
Samp3: 180 degree rotation
Samp4: 270 degree rotation
Samp5: 360 -> 0 degree rotation
Samp6: 90 degree rotation
.
.
.

For a pure 90 degree phase shift all you need to do is:

Ir = -Qi,
Qr = Ii

at every sample.

It might be helpful for you to read up on cordic rotation. Although you might not be interested in the cordic algorithm it will provide a good overview of rotations and frequency shifts.


Advantages and disadvantages of Hilbert-Huang Transform versus Fast Fourier Transform? Which one is better for frequency detection? In what circumstances? Advantages and disadvantages of both? Which would be better for a Voice Activation Detection algorithm? Which would be better for an instrument tuner? How would one detect pitch user an HHT transform? Is it faster than FFT? How would one implement HHT?

Fourier algorithm is actually a global transform that can not reflect the damping and local specialty. Unfortunately, these features are fundamental and important for characterize a nonstationary signal. on the other hand in the Fourier algorithm the estimation of frequencies are sensitive to noise. so HHT is better in your case if you are going to deal with noisy signals.


OFDM coarse freq estimation I am trying to read OFDM synchronization paper by Schmidl and Cox. The math is slightly overwhelming.

What is the formula that they use in implementation that distinguishes between coarse estimation and fine estimation?

Secondly once fine frequency offset is done, how do you determine that you need to do coarse frequency synchronization or not?  What is it that you will check to determine it? If anybody could explain with a set of equations that would  be very helpful.

The key part of the paper seems to be section IV.  There, they estimate the frequency offset by:

where $P(d)$ is given by:

And then $\widehat{\Delta f}$ is given by

if $|\hat{\phi}| < \pi$, or a correction is added:

otherwise.


How to solve this entropy integral? I am having the entropy integral below where $\mathbf x$ is a $N$ dimensional Gaussian vector having variance as $\mathbf P$ and mean zero,

$$ -\int\frac {\exp\left(\frac{1}{2}\mathbf x^T\mathbf P^{-1}\mathbf x\right)}{{(2\pi)}^{\frac{1}{2}}|\mathbf P|^{\frac{1}{2}}}\log_2\left(\frac {\exp\left(\frac{1}{2}\mathbf x^T \mathbf P^{-1}\mathbf x\right)}{{(2\pi)}^{\frac{1}{2}}| \mathbf P|^{\frac{1}{2}}}\right) d\mathbf x$$

How do I solve this integral ? My approach is,  let

$$
\mathbf x=\begin{bmatrix}
    x_{1} \\
    x_{2}  \\
    \vdots\\
    x_{n}
\end{bmatrix}\quad\text{and}\quad \mathbf P=\begin{bmatrix}
    \sigma_{1} & 0 & 0 & \dots  & 0 \\
    0 & \sigma_{2} & 0 & \dots  & 0 \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & 0 & \dots  & \sigma_{n}
\end{bmatrix}$$

I am getting an integral like
$$-\displaystyle\int \frac {\exp\left(\frac{1}{2}\left(\frac{x_1^2}{\sigma_1^2}+\frac{x_2^2}{\sigma_2^2}+\ldots+\frac{x_n^2}{\sigma_n^2}\right)\right)}{{(2\pi)}^{\frac{1}{2}}|\mathbf P|^{\frac{1}{2}}}\log_2\left(\frac {\exp\left(\frac{1}{2}\left(\frac{x_1^2}{\sigma_1^2}+\frac{x_2^2}{\sigma_2^2}+\ldots +\frac{x_n^2}{\sigma_n^2}\right)\right)}{{(2\pi)}^{\frac{1}{2}}|\mathbf P|^{\frac{1}{2}}}\right) d\mathbf x$$
but I am confused how to do this $d\mathbf x$ stuff when $\mathbf x$ is a vector. Am I doing it right?

The answer given is

  $$\frac{n}{2}\log_2(2 \pi e)+ \frac {1}{2}\log_2|\mathbf P|$$

The integration over a vector $\mathbf x$ as given by you means something like

$$\int f(\mathbf x) d\mathbf x=\int\dots\int f(x_1,\dots,x_n)dx_1\dots dx_n$$

So, you integrate over each component separately. Now, to solve your integral, you should first split the $\log$ into a difference and evaluate both integrals separately (exploiting the fact that $\exp(ab)=\exp(a)\exp(b)$. Then you will get a lot of similarly structured integrals that depend only on $\sigma_i$. If you solve them, you'll get the solution.


Estimating the bandwidth of various signals I am trying to estimate the bandwidth of different signals, with different modulations, as seen in a PSD graph.  I've heard about a "$3\ \rm dB$ rule" in this context.  I don't need extreme accuracy.  Might drawing the edge of the channels at $3\ \rm dB$ down from the top work for me?

Another way would be fitting each of your signals' PSD with a Gaussian probability density function whose standard deviation $\sigma$ and mean $\mu$ you estimate, and use the full width at half maximum (FWHM) of this approximation as your half-power bandwidth. This is known accurately with the following closed form:
$$
\mathrm{FWHM}=2\sqrt{2\ln(2)} \sigma\approx 2.355\sigma
$$
This approach is often used in practice to measure things like bandwidth expansion.


How does sampling rate of $x[n]$ relate to sampling rate in frequency domain after DFT transformation? I've got an analog signal $x(t)$ sampled at frequency $F_s$ to obtain samples:

$$ x[n] = x(t) \bigg|_{t=n/F_s} $$

I transform this signal with DFT defined as:

$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-i2\pi kn/N}$$

Is there such a thing as a sampling frequency / rate in the frequency domain? How does it relate to $F_s$?

A better way to think about it is the DFT is agnostic about the sampling rate.  What matters is the number of sample points.  If you look at the DFT definition, $F_s$ does not appear in it, but $N$ does twice.

The sampling rate can be used to figure out the duration of the sampling frame given the number of samples, or used to figure out how many sample points for a given duration.  Once the DFT is taken on the sample, the bin number corresponds to the number of cycles per frame.  The number of cycles per frame can be divided by the number of seconds per frame (the frame duration), or whatever unit of time is being used, to get the number of cycles per second (or time unit).

I think the equation the r b-j gave is better expressed as:

$$ f = k \cdot \frac{F_s}{N} $$

Where the units are:

$$ \frac{cycles}{timeunit} = \frac{cycles}{frame} \cdot \frac{\frac{samples}{timeunit}}{\frac{samples}{frame}} $$

The latter factor can be considered a conversion factor for the frequency.

When $ k = N/2 $:

$$ f = \frac{N}{2} \cdot \frac{F_s}{N} = \frac{F_s}{2} $$

Which confirms r b-j's statement about the halfway bin corresponding to the Nyquist frequency.

Hope this helps.

Ced


Does SWT/ISWT require intermediate approximation coefficients to represent/reconstruct the original signal? Taking into account that Stationary Wavelet Transform (Algorithme à trous) is not an orthogonal transform do we need intermediate approximation coefficients for signal decomposition?

For example, in maximally decimated DWT we have the following decomposition tree:

And according to the conclusions from the MRA (multi resolution analysis) to fully represent (and reconstruct) the signal it is suffice to have only the last approximation subband plus all of the difference subbands.

Does this also hold true for SWT?

If I understand well, no, you don't need intermediate approximation coefficients, as long as you have (not) sampled properly. So yes, having the  last approximation subband plus all of the difference subbands suffices.

Moreover, there are several inverses, since the transform is redundant, which can reveal better than the obvious one. Note that a stationary wavelet packet transformation  would require intermediate approximation coefficients.

If I may, orthogonality is not the question here: you could perform SWT with non-orthogonal wavelet filters. With insight from your other question Additive white gaussian noise and undecimated DWT, I propose to recall that:

a filter-bank is called critical when the number of samples and coefficients is the same (outside border). It does not need to be orthogonal
such a linear transform can generally be implemented with a square matrix
in some cases, the matrix is orthogonal, a subset of critically-sampled transforms.

In the literature, one sometimes finds "orthogonal" as a proxy for "square" & orthogonal, which is correct, but "square" being more important in that case,  as compared to redundant schemes.


How to use ImageNet Database in Matlab How to use ImageNet Database in Matlab for object classification or recognition? I didn't find any tutorial which could practically tell me how to train images on ImageNet database and classify the objects.

Use MatConvNet. They have samples for image-net classification, such as this. Pre-trained models already exist. You could also train on your own.

Check cnn_imagenet function in MatConvNet. You will find it useful.
There are tutorials. For example this one.


Handling a number of simultaneous measurements in Extended Kalman Filter? If I have a number of sensors whose measurements arrive at the same time - how can I handle them properly in the Kalman filter?

P.S.
The measurements are not necessarily taken at the same time but I only receive them batched together.

If your measurements are batch, so is your filter.

You can do a number of things.

You can sequence your measurements in natural order and update the filter for each.

You can do the previous and then approach the problem as a fixed lag smoother.

There is a literature on out of sequence and delayed updating as well, like:

  Ensemble Kalman filter for multisensor fusion with multistep delayed
  measurements Sirichai Pornsarayouth; Masaki Yamakita

  2011 Aerospace Conference Year: 2011 Pages: 1 - 10 Cited by: Papers
  (1) IEEE Conferences


Should scipy.signal.coherence be 1 for single input and output signals? I am trying to calculate the coherence between input and output signals. I thought I could work with a single input and a single output time series and calculate the coherence $\gamma^2$ between them. I did so with the scipy.signal.coherence function (See code and results below).
Now, a colleague challanged my procedure, saying that the coherence needs to be calculated by averaging several input and output signals, since otherwise the coherence will simply be $\gamma^2=1$. I wanted to prove otherwise, after all I used a scipy function and I regard those as trustworthy.
So I started out with the general definition of the coherence as [1][2] [3]:
$$
    \gamma^2 = \frac{|\mathbb{E}[X(-j\omega)Y(j\omega)]|^2}{\mathbb{E}[|X(j\omega)|^2]~~\mathbb{E}[|Y(j\omega)|^2]}, \quad \hat{\gamma}^2 = \frac{|\overline{S}_{XY}|^2}{\overline{S}_{XX}\overline{S}_{YY}},
$$
where $\overline{\cdot}$ indicates a mean value, $|\cdot|$ the magnitude of a complex number, ${S}_{XX}$ and ${S}_{YY}$ are the energy spectral density for $X$ and $Y$ respectively and ${S}_{XY}$ is the cross spectral denisty between $X$ and $Y$. And $X(j\omega)$ being the Fourier Transform of a time series $X(t)$ and $X(-j\omega)$ being its complex conjugate. Let's define those as:

$$
X(j\omega) = x_1+ix_2, \quad X(-j\omega) = x_1 -ix_2, \quad Y(j\omega) = y_1+iy_2.
$$

Now, assuming I do not have several signals, but only a single input signal $X$ and single output signal $Y$, I drop the mean $\overline{\cdot}$ and end up with the following:
$$
    \gamma^2 = \frac{|{S}_{XY}|^2}{{S}_{XX}{S}_{YY}} = \frac{|X(-j\omega)Y(j\omega)|^2}{|X(j\omega)|^2|Y(j\omega)|^2} = \frac{\left(|X(-j\omega)||Y(j\omega)|\right)^2}{|X(j\omega)|^2|Y(j\omega)|^2}= \frac{\left(\sqrt{x^2_1+x^2_2}\sqrt{y^2_1+y^2_2}\right)^2}{\left(x^2_1+x^2_2\right)\left(y^2_1+y^2_2\right)}=1.
$$

In case I did not make a mistake, my scipy.signal.coherence function should give me a coherence = 1. But as shown below, it does not. Instead the coherence is close to zero everywhere, except at the signals' natural frequencies. Please note that I get $\gamma^2=1$ only if I chose $X=Y$.

I must be missing something here but I just don't see it. Any ideas?

Many thanks in advance!
Max

Result:

Code:

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
from scipy.fftpack import fft

fs = 5e3 # sampling frequency
N = 10000 # number of observations
T = N/fs # max time T
amp = 20 # sine wave amplitude input
amp2 = 2*amp # sine wave amplitude output
freq = 500.0 # sine wave frequency input
freq2 = 500 # sine wave frequency output
time = np.arange(N) / fs # time vector t_0 up to T
# create data
x = amp*np.sin(2*np.pi*freq*time) # input
y = amp2*np.sin(2*np.pi*freq2*time)+np.random.normal(scale = 10, size = len(x)) + 100 # output

fig, axes = plt.subplots(1,3, figsize=(15,5))

# plot time signals
axes[0].plot(x[:400], label = 'Input', alpha = 0.7)
axes[0].plot(y[:400], label = 'Output', alpha = 0.7)
axes[0].legend()
axes[0].grid(linestyle='--')
axes[0].set_xlabel('Time')

# plot fft input
xfft = fft(x)
N = len(x) # number of sample points
T = 1 / fs # sample spacing
xfreq = np.linspace(0.0, 1.0/(2.0*T), N//2)
axes[1].plot(xfreq, 2.0/N * np.abs(xfft[0:N//2]), label = 'Input', alpha = 0.7)

# plot fft output
yfft = fft(y)
N = len(y) # number of sample points
T = 1 / fs # sample spacing
yfreq = np.linspace(0.0, 1.0/(2.0*T), N//2)
axes[1].plot(xfreq, 2.0/N * np.abs(yfft[0:N//2]), label ='Output', alpha = 0.7)
axes[1].legend()
axes[1].set_xlim(1,1000)
axes[1].grid(linestyle='--')
axes[1].set_xlabel('Frequency [Hz]')

# plot coherence:
f, Cxy = signal.coherence(x, y, fs)
axes[2].semilogy(f, Cxy, color='black')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Coherence')
axes[2].grid(linestyle='--')

plt.show()

Your colleague is correct.

In

  Carter, G. Clifford. "Coherence and time delay estimation."
  Proceedings of the IEEE 75.2 (1987): 236-255.

a closed form expression for the pdf of the estimate $\mid \hat{\gamma^2}\mid$, the  $p(\mid \hat{\gamma^2}\mid |\; N,\mid \gamma^2 \mid)$, is given and the confidence pdf for $N=1$ is flat on $[0,1]$


Scalogram (and related nomenclatures) for DWT? My understanding of the scalogram is that, for a particular row, the scores of the projection of the input signal with the wavelet at a particular displacement is shown. Across rows, the same thing applies, but for dilated version of the wavelet. I thought that scalograms can be defined for all types of wavelet transforms, that is, for the:

Continuous wavelet transform
Discrete wavelet transform
Redundant wavelet transform

However upon further investigation is seems that the scalogram is only definable for the CWT. Based on this I have multiple inter-related questions that google has not sufficed for ATM.

Questions:

Is it true that the scalogram is not defined for the DWT or RWT? If so, why not?
Let us say an $N$ length signal has a 10-level decomposition by using DWT. If all levels are plotted as an image, (that is, a $10xN$ image), what is this image called?

As an example of a DWT 'scalogram', here is one for AWGN:

Concerning the same signal, suppose we instead plot the approximation MRA of the signal at all levels. (So again, a $10xN$) image. What is this image called in proper terminology?
For example, here I have shown approximation MRAs and detail MRAs for AWGN. (Clearly they are not the same as 'scalogram' of DWT).

Thanks!

Continuous wavelet transform is suitable for a scalogram because the analysis window can be sized and placed at any position.  This flexibility allows for the generation of a smooth image in both the time in scale (analogous to frequency) directions.  The continuous wavelet transform is a redundant transform because the analysis window can overlap.  In fact the CWT is considered infinitely redundant.
Discrete wavelet transform is a non-redundant transform.  It was developed so there would be a one to one correspondence between the information in the signal domain and the transform domain.  This tight correspondence makes the DWT more suitable for use in signal reconstruction.  The analysis windows are fixed in both the time and scale directions, so if you plot the resulting DWT coefficients you will end up with a grid of boxes that start out large at one end of the scale axis and end up small at the other end.  This representation isn't very satisfying for visual analysis of a signal.  It certainly can be done, but I haven't seen anyone bother to do it.  The plot is also referred to as a scalogram.
Redundant Wavelet Transform:  I had no previous experience with this, but thanks to comments from the OP, I found that the RWT or Stationary Wavelet Transform (SWT) is a discrete wavelet transform that has redundancy introduced to make the transform translation invariant.  Furthermore, I found a reference that does a nice comparison of transform types as they apply to speech analysis.  In this article, the transform results are all plotted and for any case of wavelet transform, the plots are all referred to as scalograms (this includes the DWT and a version of RWT).  You can see how the various transform types present themselfs visually in the article.  For reference, here is a link to the article:  http://www.math.purdue.edu/~lipeijun/paper/2005/End_Gen_Li_Fra_Sch_JASA_2005.pdf

MRA - My encounter with this term is in association with multiresolution analysis.  This applies to all wavelet transform types, but usually is discussed in the context of the DWT and its realization as a set of filter banks.  In this context the result of a MRA is the same as the result of a DWT and the plot of such results (a plot of a set of numbers) would still be a scalogram.  Here is another paper that discusses MRA: http://alexandria.tue.nl/repository/books/612762.pdf

The following is an example of CWT and DFT Scalograms:


How to design a FIR low pass filter and force a subset of filter coefficients to 0? I'm looking into interpolating a signal by 2x and came into this matlab function interp:

  interp

  Interpolation — increase sampling rate by integer factor

  Algorithms

  interp uses the lowpass interpolation algorithm 8.1 described in 1.

  1.It expands the input vector to the correct length by inserting 0s between the original data values.

  2.It designs a special symmetric FIR filter that allows the original data to pass through unchanged and interpolates to minimize the
  mean-square error between the interpolated points and their ideal
  values. The filter used by interp is the same as the filter returned
  by intfilt.

  3.It applies the filter to the expanded input vector to produce the output.

Then about intfilt:

  intfilt

  Interpolation FIR filter design

  Algorithms

  The bandlimited method uses firls to design an interpolation FIR
  filter. The polynomial method uses Lagrange's polynomial interpolation
  formula on equally spaced samples to construct the appropriate filter.

Finally the documentation of firls is here

In short, the algorithm of interp is to insert l-1 zeros between each pair of neighboring input samples (up-sampling) and FIR filter it, where l is the multiplier of interpolation. The FIR filter is designed by intfilt.

E.g. intfilt(2,4,0.5) returns a 15-tap FIR filter with 0.5*fs/2 cut-off for 2x interpolation, in this case:
    [-0.0068, 0.0000, 0.0395, -0.0000, -0.1427, 0.0000, 0.6098, 1.0000, 0.6098, 0.0000, -0.1427, -0.0000, 0.0395, 0.0000,-0.0068]

Observation:

In the above example: all the even coefficients are 0 and the middle coefficient is 1. Recall that this FIR filter is to be used with a 2x up-sampled signal, this means the "original" sample in the input will be unchanged (only multiplied with the center 1), while all the inserted 0s will the calculated by the FIR filter.

This offers two useful properties: 1) computational efficiency: "original" samples are unchanged; 2) lossless-ness: interpolation then down-sample leaves the original signal unchanged.

Question:

I found the source code for intfilt and it calls h=firls(n-1,F*2,M); with n=15, F = [0, 0.1250, 0.3750 , 0.5000] and M=[2, 2, 0, 0], then returned h is the same as h = intfilt(2,4,0.5);

Why this n/F/M combination would set all even coefficients of h to 0?

The filter you describe is called a halfband filter. It is symmetric with respect to half the Nyquist frequency ($f_s/4$). The function firls computes the least squares optimal filter. As long as the desired magnitude and the edge frequencies result in a symmetric specification, the optimal approximation will also be symmetric, and, consequently, will be a halfband filter with every other sample away from the center tap equal to zero (up to numerical accuracy).

So any specification for firls of the form

f = [0,0.5-x,0.5+x,1];
m = c*[1,1,0,0];

with $0<x<0.5$ and $c>0$, and an odd filter length will result in a halfband filter.


image feature extraction We are given a set of grayscale image patches obtained from images after edge detection. Each patch is 10x10 pixels with intensity varying between 0 and 255 for each pixel. This set may contain a very few (maybe one or two) or a large number of patches. We want to represent the set by a single 10x10 pixel patch or “feature” such that the feature captures the invariance or common properties of the set. Intuitively speaking, the feature should be the intersection of the edges in the image patches in the set. For example, if all the patches contain a vertical edge at a fixed position and also contain other edges in different orientations and positions, then our feature should contain only the vertical edge at the fixed position and no other edge. Unfortunately, the operation of intersection is not well-defined for images. How can we approach to compute such a feature for a set of image patches?

I may be wrong if i have not understood the question! I am trying to give a rather elementary introduction here. I can refine things and be more rigorous as suited.

What you are looking for is that of 100 (or 1000) patches, which patch is the most representative patch of all.

For simplicity if the size of a patch is 1x1. So it is just a scalar. In this case, you can just find a simple average $RepScalar = Avg = \sum P_i$.

You can also apply median rather than averaging if appropriate. A more generalized way is to treat this as random variable and find expected value as $E(P) = \sum PDF(x) \cdot x $

I don't know if you are really conversent with this: so i am leaving you with a tutorial on this: http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/Chapter6.pdf

Now in your case the $P$ is a vector instead of scalar.

So the averaging will become

$$
RepPatch[i] = \sum_{k} P_k[i]
$$

where $k$ is a patch number and $i \in { 1, ... 100 }$ for $10x10$ patch size.

If you know the PDF or if you can estimate PDF you can find Expected value of this but for a vector of size of 100 (10x10) it is too complex.


How to implement the imaginary part of FFT equation applied to a signal in computer/C++? I am trying to apply FFT on a time domain signal in C++, however I am having trouble with the imaginary part i.e square root of -1. In the FFT equation we have,

$$X =\frac{1}{2\pi}\int_a^b x(t) \ e^{-j\omega t} \ dt$$

we know that:
$$
e^{-j\omega t} = \cos(\omega t) - j\sin(\omega t) \tag{2}
$$

we also know that:

$$ \ j = \sqrt{-1}\ $$

Problem :  I am trying to apply the fft equation on a signal and the problem of implementation arises when coding the equation(2) in C++. As, we know the value of $ \ j = \sqrt{-1}\ $ is infinite or if tried to solve using a calculator, we get a Math Error*. And so is my program giving an error and I'd like to get suggestions and advises on how to solve this issue and apply the FFT equation in C++.

In computer programs, complex numbers are usually represented by an array, field or vector containing two values: one value represents the real part, the other represents the imaginary part of the complex number. The imaginary unit itself is usually not represented at all, because this information would be redundant.

An alternative representation is to store the magnitude and phase angle of a complex number.

Specifically in C++, there is the class complex<T> from the STL library that takes care of storing two values of type T for every complex value. You should use that, it's part of the standard libraries.


Morlet wavelet time and frequency resolution I want to calculate Morlet time and frequency resolution. the Morlet wavelet function is define as : $\psi(t)=\frac{1}{\sqrt{\pi f_b}}e^{j2\pi f_c}e^{-t^2/f_b}$
Note: I know the answers, but I don't know how to achieve it on my own.I want someone to explain the time and frequency resolution relation to me please. here are some hint that I got from paper and might help you.
in some paper I found that Morlet wavelet time and frequency resolution is :
 $\Delta t=\frac{f_c\sqrt{f_b}}{2f_i}$
 $\Delta f=\frac{1}{2\pi f_c \sqrt{f_b}}$
Because generally, the wavelet time and frequency resolution is define as:
$\Delta t=s \Delta t_\psi$   and    $\Delta f= \frac {\Delta f_\psi}{s}$  (1)
where 's' is scale. it is also evident that the frequency and scale are related to each other by $f_i=\frac{fc}{s_i}$ .   $\Delta t$ and $\Delta f$ are wavelet time and frequency resolution,respectively. $\Delta t_\psi$ and $\Delta f_\psi$ are morlet function time and frequency resolution which are as follow:
$\Delta t_\psi=\frac{\sqrt{f_b}}{2}$ and $\Delta f_\psi=\frac{1}{2\pi \sqrt{f_b}}$ (2)
the Heisenberg uncertainty principle also says: $\Delta t \Delta f >= \frac{1}{4\pi}$
Basically i want to prove (1) and (2). other things can be obtained by substituting.
Thanks

Consider shifted and scaled versions of a mother wavelet $\psi(t)$:

$$\psi_{a,b}(t)=\frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right),\quad a>0,\;b\in\mathbb{R}\tag{1}$$

By the definition of the Fourier transform

$$\Psi(\omega)=\int_{-\infty}^{\infty}\psi(t)e^{-i\omega t}dt$$

it can be shown that the Fourier transform of $\psi_{a,b}(t)$ is

$$\Psi_{a,b}(\omega)=\sqrt{a}e^{-ib\omega}\Psi(a\omega)\tag{2}$$

Note that the scaling factor $a$ appears in the denominator of the argument of $\psi_{a,b}(t)$ and in the "numerator" of the argument of $\Psi_{a,b}(\omega)$. This implies that stretching in the time domain ($a>1$) implies compression in the frequency domain and vice versa. Consequently, if $\Delta t$ and $\Delta f$ denote the "spread" of $\psi(t)$ in time and frequency, it follows from (1) and (2) that the corresponding spreads of $\psi_{a,b}(t)$ - call them $\Delta t_{a}$ and $\Delta f_{a}$ - are given by

$$\Delta t_a=a\Delta t\quad\textrm{and}\quad\Delta f_{a}=\frac{\Delta f}{a}$$

As for the definition of $\Delta t$ and $\Delta f$, there is not way to "prove" the results that you stated. It is just a matter of defining the spread of a function. Since $|\psi(t)|$ is a Gaussian with standard deviation $\sigma=\sqrt{f_b/2}$ it seems natural to define $\Delta t$ as being proportional to the standard deviation, but there is no "correct" way to choose the proportionality constant. Once $\Delta t$ has been defined, the corresponding definition of $\Delta f$ can be obtained from the uncertainty principle.


What Is the Correct Way to Apply Richarson Lucy Deconvolution to Luminance Data? My question concerns the Richardson–Lucy deconvolution algorithm, which is described in Richardson's original paper. I am interested in applying it in the context of a raw image converter for digital photography. The RawTherapee discussion that motivated this post can be found here.

By using a Poisson noise model, the original paper assumes the algorithm is being applied to the linear data recovered from the sensor, which in particular is not gamma corrected any does not have any additional nonlinear tone mapping applied. A Poisson random variable is no longer Poisson after a nonlinear transformation, so the model no longer applies to such data.

My understanding is that if we simply demosaic the data from the Bayer array on the sensor, then apply deconvolution to the resulting RGB data, we are in good shape because it is linear.

However, it seems traditional to apply the deconvolution only to the luminance data. (I'm not totally sure why, or whether it is advisable, so please consider this a sub-question.) RawTherapee does this, for example, by converting the image to a LAB color space and deconvolving the L channel.

But there's the issue. It doesn't seem like L is linear in the number of photons captured. I found some data that suggests reflectance percentages of (12.5, 25, 50, 100) correspond to L values of (42.4, 57.0, 75.8, 100). Possibly this conversion depends on the illuminant being used; I am not well-versed enough in color science to say. However, the point is that the mapping is definitely non-linear.

So, is there a practical way to apply the Richardson–Lucy algorithm to a "linearized" version of the luminance data? Is this even desirable? I am not sure what the best practices are here, so any references you could provide would be greatly appreciated. I could not find any mention of this issue in the references I looked at.

I think the answer is in Jim Kasson's Post (The Zone System and Digital Cameras) you linked to:

If you are only moderately familiar with the Zone System, you may be asking yourself about the utility of recording light brighter than 100% reflectance. The 100% reflectance calibration assumes a matte subject with perfectly Lambertian reflectance. The higher light levels occur when the subject has some specularity or gloss to it. Water is a classic, and photographically-important example. Metals, glass, plastics, skin, stone and many paints are also often far from Lambertian. So we often need the headroom to accommodate those materials.
All of the above assumes that the camera is calibrated using UniWB, and that the UniWB hack works perfectly to predict raw file values. However, UniWB never achieved much penetration among digital camera users, and I think the number of users is waning today as cameras get more and more dynamic range, so it’s worthwhile looking at what happens with widely-used camera setups.

So the model to of Reflectance and the actual number of photons the sensors is exposed to are linear only if the camera is behaving according to UniWB Model (See also Guillermo Luijk - UniWB Tutorial) and the material used is indeed perfectly Lambertian.
Both probably doesn't hold and create some non linear effects.
CMOS / CCD Sensors are indeed linear in their response (At least before they are clipped).
So, given a RAW data with Linear Demosaicing (Pay attention that some algorithms may be non linear in this step) you can indeed apply the Lucy Richardson Convolution on the data.
Pay attention that LAB Color Space is a non linear transformation of the RGB colors (even without the Gamma Correction).
So I'd apply the deconvolution in the following order:

Demosaic RAW Data.
Convert data into YCbCr Color Space (Or similar) which decompose the colors into luminosity and color data with Linear Transformation.
Apply the deconvolution on the Y channel.
Convert back to RGB based color space.


Lightweight test for bimodal distributions? Here, and in the stats stackexchange, seem to be answers that reference tests for bimodal distributions that involve iterative binning or iterative curve fitting methods.  However "eyeballing" a plot of a data set often shows a clear bimodality (say a 10 dB dip or several standard deviations  between two clear mode peaks, etc.), versus a single "hump", or something ambiguous (less than a 3 dB dip).

Are there any lightweight algorithms (computationally efficient, single pass, or deterministic low iteration count tests) which can separate out the easy cases (clearly unimodal, or clearly bimodal), versus data sets that may require more computationally intensive tests for likelihood of multi or bi-modality?

Or, in signal processing terms, how can one non-graphically and quickly test whether some area in a noisy spectrogram is more likely to contain 2-tone FSK (with an unknown separation), rather than noise or some other kind of modulated signal?

An easy way, especially one that is used in the detection and SNR estimation of PSK signals in a noisy signal, is based on stochastic moments of the received amplitude:

The kurtosis can be used as measure of how little "centered" a symmetric distribution is.

You could calculate a sample kurtosis frequency bin-wise in your spectrogram; those bins with high curtosis are likely to contain more than just Gaussian noise.

In a slightly more elaborate version of that algorithm optimized for detecting 2-FSK, you could use that to identify "candidate bins", subtract their mean power to have zero mean variables, and then calculate the covariance of pairs of these (i.e. calculate a covariance matrix). You're looking for a strong negative covariance – whenever one bin has above-average power, the other has low power.


Can you recommend an efficient encoding algorithm for sound with long periods of silence? As a side personal project, I'm coding an Android application that I plan to use to record the sounds of our first baby. I'm kind of new to the world of Android development and sound processing, so please bear with me if this question is naive :)

What I'm trying to do is to constantly save the sounds in our baby room, splitting the recordings in 2-hour intervals. I plan to keep this data for future analysis (e.g. plotting crying times, which I'll correlate with other events, and maybe cry identification).

So far the application is working fine. However, I wonder if there is an efficient codec I can use considering that 90% of the time there will be silence. I've tried modifying the PCM raw data before enters the codec setting to zero any values below a specific threshold, but this did not reduce the file size substantially. This was a wild guess, just in case this would "help" the codec (I've tried AAC/HE and MP3, not yet AMR).

Can you please point me in the right direction?...

I don't mind using a solution that requires post processing of the files. For example, I've looked into http://sox.sourceforge.net, but that strips silence periods reducing the recording length (I want to keep the timestamps and length untouched).

Note: using SoX to split the file in parts using a naming scheme that allows me to reconstitute original times it's an option... but I could not come with that script, yet).

Bonus question: will AMR-NB be enough to encode a baby cry, if I want to run some processing later? (for example this or this) or... should I stick to AMR-WB?

Thank you,

Update July 22nd:

Thanks everyone for the responses so far. The idea of using just WAV injecting zeros, and later using bzip2, etc. to compress the repetitive data is really good (thanks @MBaz).

@Fat32: I've erased the files, but I recall sizes around ~20Mb/hr, with no substantial change after replacing low PCM values (approximately 50% of the stream) with zeros.

Another update: AMR-WB is giving me better results than AAC/HE, ~6 Mb/hr. So I'll probably stick to it, however I'm not sure how the low bitrate + sampling rate will affect things if later I decide to run some infant identification analysis later (i.e. if the loss of information will remove data points crucial for the neural net).

Thresholding should reduce the file size to the limits of the encoder being used. AAC/H.264 are using the state of the art general purpose audio compression codecs (which could only be beaten by application specific optimized codecs, such as speech vocoders etc).

May be your files are already minimal? Put some numerical data such as the sampling rate, average bits per sample, recording duration and file sizes...


Techniques for detecting small blobs in noisy images I am trying to write a program that uses computer vision techniques to detect (and track) tiny blobs in a stream of very noisy images. My data is of two types: one set of images are not so noisy, which I am just using to try different techniques with, and the other set are noisier, and this is where the detection needs to work at the end. Here are some samples of the higher-noise images. The blobs that need to be detected are the small black spots near the center of the image.

Initially I started off with a simple contour/blob detection techniques in OpenCV, which were not very helpful. Eventually I moved on to techniques such as "opening" the image using morphological operators, and subsequently performing a Laplacian of Gaussian blob detection to detect areas of interest. This gave me better results for the low-noise versions of the images, but fails when it comes to the high-noise ones: gives me too many false positives. Here is a result from a low-noise image (please note input image was inverted).

The code for my current LoG based approach in MATLAB goes as below:

while ~isDone(videoReader)
    frame = step(videoReader);
    roi_frame = imcrop(frame, [660 410 120 110]);

    I_roi = rgb2gray(roi_frame);
    I_roi = imcomplement(I_roi);
    I_roi = wiener2(I_roi, [5 5]);
    background = imopen(I_roi,strel('disk',3));

    I2 = imadjust(I_roi - background);
    K = imgaussfilt(I2, 5);
    level = graythresh(K);
    bw = im2bw(I2);

    sigma = 3;
    % Filter image with LoG
    I = double(bw);
    h = fspecial('log',sigma*30,sigma);
    Ifilt = -imfilter(I,h);

    % Threshold for points of interest
    Ifilt(Ifilt < 0.001) = 0;
    % Dilate to obtain local maxima
    Idil = imdilate(Ifilt,strel('disk',50));

    % This is the final image
    P = (Ifilt == Idil) .* Ifilt;

Is there any way I can improve my current detection technique to make it work for images with a lot of background noise? Or are there techniques better suited for images like this?

I'd recommend looking at the state of the art using deep convolutional neural networks (CNNs) for object localization and detection.

The difference is that a deep CNN architecture facilitates a data-driven approach, where the output of your CNN is a trainable classifier, that with more and different kinds of noisy/scaled/translated/permuted input data, can learn the feature maps that allow the CNN to generalize much better (say to different light & noise conditions) than your hand-engineered features (morph filtering, laplace of gaussians, thresholding) used with a static threshold (0.001)

Start on github.


Why do poles in the left half of the S plane make a system stable? A point on the S-plane (where $s=\sigma+j\omega$) represents a signal with a given frequency (given by the imaginary component) and which either decays, increases or stays stable (depending on the value on the real component).

Doing the maths, by converting $e^{-st}$ to a cosine and sine pair, I can see that points the left hand side of the plane describe signals that increase infinitely and points in the right hand side of the plane describe points that decrease to infinity.

Given this, why is it true that having poles of a transfer function (the frequency values which make the gain of the system infinite), lying in the left hand side of the S-plane (the side which makes signals increase infinity), makes a system stable ?

An $n$th order linear system is asymptotically stable only if all of the components in the homogeneous response from a finite set of initial conditions decay to zero as time increases, or mathematically written:
$\displaystyle{\lim_{t \to \infty}}\sum_{i=1}^{n}C_ie^{p_it}=0$ ($p_i$ are the system poles).  As time increases (in a stable system) all components of the homogeneous response must decay to zero.If any pole has a positive real part there is a component in the output that increases without bound, causing the system to be unstable. So, in order for a linear system to be stable, all of its poles must have negative real parts (they must all lie within the left-half of the s-plane).  An "unstable" pole, lying in the right half of the s-plane, generates a component in the system homogeneous response that increases without bound from any finite initial conditions.


How to simulate hyper-spectral images Hyper-spectral images are images with many bands $(\ge 100)$ while the normal color images have only three bands $\rm (r,g,b)$.  It is not difficult to simulate an artificial gray-scale image with the following Matlab codes:

width = 200;
height = 200;
img = ones(height,width)*10;
img(floor(height/4):floor(height*3/4),floor(width/4):floor(width*3/4))= 200;
noise = 5+10*randn(height,width);
img = img+noise;
img(img<0)=0;
img(img>255)=255;
figure; imshow(img,[]);

The following simulation image can be obtained:

In the simulation image, the object is a square in the center and the rest is the background. The object is of gray-value 200, and both object and background are contaminated by Gaussian noise. It is also possible to perform blurring before adding noise. The blurring is used to simulate the system point spread function.

Then comes to my question: how can we simulate a hyper-spectral image?

For the time being, my solution is as follows:

Find a hyper-spectral profile from the hyperspectral database for a certain object (this kind of database exists?)
For each band get the object's gray-value from the profile, and treat this band as if it is gray-scale image when creating simulation image.

I am not sure whether this is the right way of generating simulation images. Any ideas?


Frequency Response with Delta Function? I am trying to find frequency response and magnitude of the frequency response of the following system impulse response:
$$h[n] = 2\delta [n] + 2\delta [n-1]$$

I understand, that through the DTFT:
$$H(\omega ) = 2(1 + e^{-j\omega})$$

However, I am unsure how to find the magnitude of the frequency response.  The solution is:
$$|H(\omega )| = 4 \left| \cos \left(\frac{\omega }{2} \right) \right|$$

What steps are needed in order to convert the above frequency response to the magnitude of the frequency response?

Hint: Expand $e^{j \omega}$ using Euler's formula. This will enable you to express $H(\omega) = R(\omega) + j I(\omega)$ where $R$ and $I$ are the real and imaginary parts and then the magnitude of the frequency response is just $\sqrt{ R(\omega)^2 + I(\omega)^2}$.


Dividing Audio signals into Barks Bands I want to apply the following steps in matlab

1) Read the audio data

2) Apply window function on that data

3) Apply fft on the result of window function

4) Make power spectrum of fft

5) divide power spectrum in Barks bands

[y,Fs,bits] = wavread('file1.wav');
Nsamps = length(y);
ham= hamming(Nsamps);
ham_fft = abs(fft(ham));
ham_fft = ham_fft(1:Nsamps/2);      %Discard Half of Points
plot(f, ham_fft)
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of Hamming Window of  a an audio');

I am not getting the correct waveform , want to draw the spectrum and then Divide into barks bands .

Read the audio data

[y,Fs,bits] = wavread('file1.wav');

Apply window function on that data

ham=hamming(Nsamps);
windowed = y .* ham;

Apply fft on the result of window function

ham_fft = fft(windowed);
ham_fft = ham_fft(1:Nsamps/2);      %Discard Half of Points

Make power spectrum of fft

PowSpec = abs(fft(ham_fft)).^2;

divide power spectrum in Barks bands

So here maybe start interesting things, based in wikipedia the 24 Critical Barks Bands in Hertz is:

20,100,200,300,400,510,630,770,920,1080,1270,1480,1720,2000,2320,2700,3150,3700,4400,5300 6400,7700,9500,12000,15500

Now to split the power spectrum in Barks bands you need to do something like:

BandBarks = [20 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];

for l = 1:length(BandBarks)-1
    CurrentBand = [BandBarks(l):BandBarks(l+1)];
    PowSpec(CurrentBand,:) %% here is the power spectrum for each Barks bands
end

OK by your commented error seems that your audios are in stereo, convert to mono, or get just one channel!

my example above assumes that each BIN number is the equivalent frequency (course this is nothing common was just an example), to find what bin is your frequency you need calculate:

$$ BandBarks * (\frac{Fs}{length(FFT)})$$

So based in this equation:

Bins = zeros(length(BandBarks),1);
for j = 1:length(BandBarks)

    Bins(j)=round(BandBarks(j)/(Fs/length(FFT)))

end

Now here how to plot individually each Bark band:

for l = 1:length(Bins)-1
    CurrentBand = [Bins(l):Bins(l+1)];
    figure(l);
    plot(PowSpec(CurrentBand,:));
end

I think I can not be clearer than that


Speech / Audio Separation I'm looking for a way to implement a Speech Separation from a .wav file. Let's say I have a .wav signal containing the following words:

"The cat flew away"

Now I don't really want to do the speech separation after I've split the signal inside blocks of size NxM because I feel as though I could miss parts of the word.

I have thought about the following algorithm:

1) When the wave signal is being processed, I process the signal and remove any signals that do not warrant consideration

2) When the signal starts to increase, we can therefore infer that this is a word and this get's counted. When the signal begins to decrease, stop counting as this could be white-noise.

The algorithm above would therefore infer that there is three words and then I can perform all of the actions needed on the signal..

Do you think this is an appropriate algorithm OR could anyone recommend any papers etc.. Which would enable me to do this?

EDIT:

Assume that I have a signal that contains "calls" from different bats. These signals have been recorded. What I basically need to do isolate the different bat calls from a signal. So, I would have the following:

-> Bat1 -> Bat2 -> Bat3, ....

So, it would work like speech segmentation. But, instead of segmenting for words, I would segment for Bat calls. Because, from what I believe, each "bat call" will be of a different or similar frequency. So therefore, would it be possible to some how "split the signal into blocks" and then calculate the frequency of the block, see if this matches with the frequency vector, and, if it does.. Keep it, otherwise remove it?

I can advice this ones :
http://www.cse.ohio-state.edu/pnl/software.html
There are papers and codes.

Speech seperation is not an easy issue I advice you read some books before.
This book which you can check its examples in the book is helping me much:
http://auditoryneuroscience.com/

The thing you said when the signal start can be one feature ,but Speech sounds change
constantly in a manner that is fast, but not too fast, and somewhat unpredictable yet
not entirely irregular. What I mean you will get false negatives with only this feature.

There are other things like periodicity, amplitude, pitch etc...
I think trying draw a spectrogram of speech and examining it is the best start.


Are there any filters that look (vaguely) like this? I'm trying to determine the effect of a particular non-linear process on my data.  From simulations, I can measure the power spectra of "ideal" and "processed" images and take their ratio, which (I think) should correspond to some sort of filter or filtering process.

So, are there any filters, or filtering processes, that look (vaguely) like the attached plots?  The two plots are the same data in linear-linear (left) and log-log (right) space.

More detail: I have a pipeline that takes timestream data where each point on the timestream is associated with a point in an image.  Some "cleaning" (PCA subtraction) is performed on the timestreams before creating the images.  My goal is to determine the effect the timestream processing has on the resulting image, with the end goal being to replicate the effect on images for which no comparable timestream data is available.  The data is coming from bolometer arrays and these are images of the sky at millimeter wavelengths.

Here are the simulated images.  Left is before, right is after processing.  The spiky points around the edges are where there is no timestream data.

See also this ipython notebook http://nbviewer.ipython.org/4393835/ for a little analysis.

Since it is entirely possible that the timestream filtering process cannot be represented by any filter in image space, an answer addressing "how do I approximate these filter functions"  is appropriate.

I would approximate the filter using an finite impulse design where the wrights of the FIR filter are computed based on the spectrum that the filter passes.  You can do this by choosing samples of the desired frequency response (from the graphs above) and taking the inverse discrete Fourier transform of the samples.  The result will be the impulse response of the desired filter.  The filter is then realized as a weighted sum using the impulse response sample values as the weights.  So your filter will be $A_0 + A_1Z^{-1} + A_2Z^{-2}...$  Where $A_0$ is the first value of the impulse response, $A_1$ is the second and so on.  $Z^{-n}$ represents a sample delay of n (coming from a Z domain representation of the filter structure).  This approach is sometimes referred to as the "frequency sampling" approach to FIR filter design.  I have used this approach to create time sequences that contain specific spectral content.  Google FIR filter design, frequency sampling.


Help on implementing an IIR filter I am new to signal processing and and having difficulty understanding how to implement an IIR filter.  The reference is athttp://www.dxatlas.com/rocky/Advanced.asp in the Band scope section.  It concerns applying an IIR filter to a power spectrum with this equation:

FilteredSpectrum[i] = FilteredSpectrum[i]*(1-gain) + NewSpectrum[i] * gain

I understand that this is a classic 1st order IIR filter but to implement this, are you working data point by data point or are you using one complete spectra and then waiting for the next one (NewSpectrum) and calculating a whole new spectrum from the existing one and the new one that just came in?
I tried it this way, waiting to obtain the data for a whole spectrum, then waiting for the next one to come in and then doing the calculation. I then plot the resulting spectrum.  This then becomes the FilteredSpectrum while I wait for the next NewSpectrum.
Is this being done in the correct way?  I programmed this but saw no effect on the noise in the spectrum contrary to the article.
Can anyone tell me where I am going wrong?
Thanks

The spectrum NewSpectrum appears to be the output of the FFT of the current data block, and FilteredSpectrum on the right-hand side of the equation is the filtered spectrum of the previous block. The index i is the frequency index, so you use this recursion for each frequency point. Note that the value of gain is crucial. It should be in the range $(0,1)$, where a value of $1$ means no filtering at all, and $0$ means no update at all.


Impulse response rest conditions of LTI systems Why the impulse response (initial) rest conditions are taken like $h(t)=0$ and $h'(t) =1$, if the order of differential equation is $2$.

Consider the second order LCCDE related to an LTI system with initial rest is given as:

$$a_0 y^{''}(t) + a_1 y^{'}(t) + a_2 y(t) = x(t) $$

The solution of the given equation for the input $x(t) = \delta(t)$ is defined as the impulse response of the LTI system and denoted as $h(t)$. The formal solution of the associated second order differential equation in time domain requires that a homogeneous and a particular solution be found and are added to yield the total solution; $$h(t) = y_h(t) + y_p(t)$$

Since the particular input is an impulse, then the particular solution is zero and therefore $$h(t) = y_h(t)$$ which is computed based on the initial conditions. For a second order equation you need two initial conditions and these are derived based on initial rest assumption so that:
$$y_h(0) = 0$$ and $$y_h^{'}(0) = 1/a_0$$.

Since most typically $a_0 = 1$ then $y_h^{'}(0) = 1$. This result comes from integrating the differential equation from $t=0^{-}$ to $t=0^{+}$ while recognising the fact that the output, $y_h(t)$, cannot involve any impulses in its defintion. In fact neither $y_h(t)$ nor any derivate of it upto (but not including) second order can include any impulses in their definition, because otherwise, then on the right side of the equation derivatives of impulse should have appeared, which is not the case, by definition, however. Hence

$$\int_{0^{-}}^{0^{+}} \left( a_0 y_h^{''}(t) + a_1 y_h^{'}(t) + a_2 y_h(t) \right)dt = \int_{0^{-}}^{0^{+}} \delta(t)$$

$$a_0 \int_{0^{-}}^{0^{+}} y_h^{''}(t) dt + a_1 \int_{0^{-}}^{0^{+}} y_h^{'}(t)dt + a_2 \int_{0^{-}}^{0^{+}} y_h(t)dt = 1
$$

Now because of the intuition based on the previous paragraph, the second and the third integrals turn out to be zero (since they don't include impulses, their integrals at a point from $0-$ to $0+$ iz zero) but the first integral includes an impulse and it will be nonzero as:

$$a_0 \int_{0^{-}}^{0^{+}} y^{''}(t) dt = a_0 ( y_h^{'}(0^{+}) - y_h^{'}(0^{-}) ) =  a_0 y_h^{'}(0^{+}) = 1
$$ and it follows that

$$y_h^{'}(0^{+}) =  \frac{1}{a_0} $$

since $y_h(t) = h(t)$ and since we can replace $0+$ with $0$, and taking $a_0 = 1$ as usual, then it follows that:
$$ \boxed{ h'(0) = \frac{1}{a_0} = 1 } $$

As the required initial condition on th impulse response $h(t)$.
Similar argumentation follows for systems of higher orders.


Is there a difference between filtering a signal before or after differentiating it? I have a time series and I want to apply:

a differentiation
a Butterworth filter

Does the order theoretically (mathematically) make any difference? Does it make any difference in real life when I use numpy?

Thanks in advance!

So, first of all:

A differentiator is really just a high-pass filter. Digitally (and this is dsp.SE, so I presume this is the case), that typically means a differentiator is really just FIR with taps $[1;-1]$. A butterworth filter being a low-pass IIR, the combination of both does sound like you might want to build a bandpass, or are doing something that doesn't make very much sense; see the comparison of a very relaxed second order butterworth and a differentiator magnitude responses below:

now, mathematically, that's identical:

In time domain, you convolve ($*$) the signal $x[n]$ with the filter impulse response:

$$\begin{array}\\
y[n] &= x[n] * h_{diff} * h_{butter}\\
&\text{or, in freq. domain, thanks to the convolution theorem}\\
Y[z] &= X[z] H_{diff} H_{butter}\\
&\text{multiplication is commutative}\\
&= X[z] H_{butter} H_{diff} \\
&\text{but it's also associative}\\
&= X[z](H_{butter}H_{diff})\\
&= X[z] \mathcal Z\{h_{butter}*h_{diff}\}
\end{array}$$
which implies you can even without problem combine the differentiator and the low pass filter. To little surprise, you get the discrete derivative of the butterworth filter.


What is the algorithm behind OCR? I'm curious about the algorithm and basic principles behind Optical Character Recognition.

Can anyone explain:

The most basic OCR algorithm that can be fairly easily implemented from scratch? Or the common denominator between different OCR algorithms?
Extensions/changes made to tailor the algorithm for application specific goals?

Also see my follow up question on OCR.

Probably the simplest OCR algorithm would be the case where the characters are from a known font, in perfect rows with no distortion or rotation, and then you can find occurrences of each character in the text by doing cross-correlation with the known characters from the known font.

Or even simpler if you can create the font before it's used, with cheaper ways to recognize characters, like the MICR fonts or OCR fonts.


Calculating energy spectral density function of half-cosine function I have the following function:

$$ g(t) = A\cdot \textrm{rect}\left(\frac{t}{T}\right)\cos\left(\frac{\pi  t}{T}\right)$$

How do I find the energy spectral density function? I think? it is defined as :

$$ \int\limits_{-\infty}^\infty \left|g(t)\right|^2\, dt $$

After squaring it, I'm a bit lost since I do not know how to continue...

For a "deterministic finite energy" continuous time signal $x\left(t\right)$, as you have exemplified, the energy density spectrum, or equivalently energy spectral density, which is a real nonnegative signal, meant to be used as a function showing the energy distribution of the signal with respect to the frequency is defined to be: $$S_x\left(\omega\right) = \left|X\left(j\omega\right)\right|^2 $$  where $X\left(j\omega\right)$ is the Continuous-Time Fourier Transform $\left(CTFT\right)$ of the signal $x\left(t\right)$ given as: $$X\left(j\omega\right) = \int_{-\infty}^{\infty}{x\left(t\right)e^{-j\omega t} dt}$$

Now coming to the specific sample of your question, you should first find the CTFT $G\left(j\omega\right)$ of $g\left(t\right)$: That your signal is a multiplication of two other signals you should use a fundamental property of CTFT namely the ${multiplication}$ in time equals ${convolution}$ in frequency. Hence

$$x\left(t\right)y\left(t\right) \leftrightarrow \frac{1}{2\pi} X\left(j\omega\right) \star Y\left(j\omega\right) $$

From a table of FTs (or by carrying out the integral if you can) find that for
$x\left(t\right)= \textrm {rect} \left( \frac {t}{T} \right)$ which is defined from $t=-T/2$ to $t=T/2$, we have: $$X\left(j\omega\right) = \frac{2 \sin \left({\omega \frac{T}{2} }\right)} {\omega}$$ and for the $y\left(t\right) = \cos\left( \pi \frac {t}{T}\right)$ we have $$Y\left(j\omega\right) = \pi \delta \left(\omega - \frac {\pi}{T}\right) + \pi \delta \left(\omega + \frac{\pi}{T}\right)$$

What remains now is to compute their convolution, which is easy considering the convolution with impulse given as:$X\left(j\omega\right) \star \delta\left(\omega - \omega_0\right) = X\left(j\left(\omega - \omega_0\right)\right) $

Hence we get: $$G\left(j\omega\right) = \frac{1}{2\pi} X\left(j\omega\right) \star Y\left(j\omega\right)$$

as $$\frac{\sin \left({\left(\omega- \frac {\pi}{T}\right) \frac{T}{2} }\right)} {\left(\omega- \frac {\pi}{T}\right)} + \frac{\sin \left({\left(\omega + \frac {\pi}{T}\right) \frac{T}{2} }\right)} {\left(\omega + \frac {\pi}{T}\right)}$$

which, after simplyfying, becomes:
$$G\left(j\omega\right) =  \cos\left({\omega\frac{T}{2} }\right)  \left( \frac{-1} {\left(\omega- \frac {\pi}{T}\right)} + \frac{1} {\left(\omega + \frac {\pi}{T}\right)} \right)$$

whose magnitude is $$\left|G\left(j\omega\right)\right| = \left|\cos\left(\omega \frac{T}{2}\right)\right| \cdot \left|\left( \frac{ \frac{-2\pi}{T}}{\omega^2 - \frac{\pi^2}{T^2}} \right)\right|$$

And from which we conclude for the Energy Spectrum Density to be:
$$S_x\left(\omega\right) = \left|G\left(j\omega\right)\right|^2 = A^2 \cos^2\left(\omega \frac{T}{2}\right)  \left( \frac{ \frac{2\pi}{T}}{\omega^2 - \frac{\pi^2}{T^2}} \right)^2 $$

where I have added the last linear scaler A in $g\left(t\right)$.


Logarithmic vs Nth root transformation of Image This is an image for the explanation of Logarithmic Transformation of an Image. I want to ask considering a input pixel value rk, Value of nth root of rk should be greater than corresponding value of log(rk). But this image is not accordance with that?

Can someone explain ?

Mathematicaly, inverse logarithm (or exponential) should indeed grow faster than powers, and logarithm is slower than roots... at infinity. Here, their range is limited to finite interval, so with proper scaling, you may observe such behaviors. I will come back later with numerical examples.


Is gamma correction applied in my setup I am working on a scientific application where I would like to evaluate luminance in an image, however I am unsure if what I receive from my setup is luminance (Y) or luma (Y'). I use a camera with a TOSH IK CU51 (PAL) control unit which is connected via S-Video to a Canopus ADVC-300 analog to digital converter. The digitized data are streamed to a computer and presented in YUV2 format.

Now the problem: I read into luminance and luma and realized that I do not know if gamma compression is applied somewhere on the way between the camera and the computer, so that really I am getting Y'UV2 instead of YUV2, is there e.g. a definition for a PAL signal to be gamma compressed and for an S-Video signal not to be corrected?

Edit: As there are no answers yet, I think I should be more specific: Is there any flexibility in the PAL or S-Video standard or in how a device is generally expected to work where I would have to take gamma correction into account?
A little more background: luminance is proportional to energy per spatial area. If gamma correction were applied then I'd get a somewhat distorted value. If I knew gamma, I could adjust for this.

You should photograph a known calibration target with your system, (like this one or the one that I put in the thumbnail).

This will allow you discovering the OECF of your system. You can either calculate it by yourself, given the technical data sheet or use third-party software like Imatest, ImageAnalyzer or others.


Frequency-ratio at baseband IQ-mixing I was wondering by what logic or mathematics the frequency-ration $$ \frac{f_{lo}}{f_s} $$ in the sine-functions comes into existence?

Source

that's pretty simple: in digital domain, real-world frequencies have no inherent "meaning".

That is, a period of a real-world sine of frequency 10 Hz might be 2, 200, 1233 or whatever samples long, depending on the sampling rate.

Thus, periods in DSP can only be measured in samples; for example, a sine of period 1 second has a period of 5.4 samples if the sampling rate was 5.4 S/s.

Now, logically, that also means that frequencies can only be related to the sampling rate. Hence, your $\frac{f_{LO}}{f_\text{sample}}$ (which is a unitless thing!) oscillator corresponds to an analog oscillator of frequency $f_{LO}$ (which has Hertz as a unit).

With your question in mind, it strikes me as odd that you're looking at a figure describing a receiver that mixes in the digital domain – in software defined radio, most receivers try to minimize the rate at which the ADC needs to run, simply for cost, efficiency, complexity, ADC dynamic range, and noise bandwidth reasons, and thus have the IQ mixer on the analog side – I'd call your design a bit of a "speciality". So, if this is you learning about direct conversion receivers, you might have picked an unlucky reference. If this is you learning about different types of digital receivers: That is extremely close to the theoretical perfect Software Radio, that does all the signal processing in digital domain.

Also note, in the figure, $f_s \overset !\ge 2\cdot f_{max}$ of the low-pass filter, not only one time that frequency (I assume "GHz" just means "order of magnitude is Gigahertzes"). That's why direct sampling is really an expensive technology for high-frequency signals (and why undersampling with filters or an IF makes a lot more sense, most of the time, since you don't care about everything happening between 0 Hz and the maximum frequency of your signal of interest, usually, but only about some limited bandwidth around a center frequency).


input output read/write frequency I am trying to identify a system based on input/output response and thereby estimate a transfer function.
I generated a frequency sweep function in Mathematica which gives me the discrete values of a frequency sweep.

g[t_] := 0.799760*Sin[2 Pi (3 t/333.3 + 1) t];
Plot[g[t], {t, 0, 1000000000000000000000000000000}]
g[#] & /@ Range[0, 5, 0.001]

So now I have a bunch of data but at what rate should I input it into the system, and at what rate should I read it?

How should I apply Nyquist frequency rule here?

Best to base it on the highest frequency response you anticipate for the system and implement twice that frequency with the samples you input and the data you read.