Convolution of Two Kernels [OpenCV] I've been completely stuck on a portion of my assignment for a few days now. After plenty of searching around, I have been unsuccessful in discovering information that leads me to the correct solution. That said, this is for OpenCV in Python, using Numpy for matrix calculations.

The problem statement:

  Construct the derivative of Gaussian kernels,  and  by convolving the above two kernels: =∗; =∗

Additionally, here's a screenshot of the problem in full form:

What I Have Attempted:

import cv2
import numpy as np

def my_Normalize(img):
if len(img.shape) != 2:
    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) #convert img to grayscale
img_float = cv2.normalize(img.astype('float'), None, 0.0, 1.0, cv2.NORM_MINMAX)
return img_float

def my_DerivativesOfGaussian(img, sigma):
#the two 3x3 Sobel kernels
#Sx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], np.float)
#Sy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], np.float)

sobelx = cv2.Sobel(img,cv2.CV_32F,1,0,ksize=3)
sobely = cv2.Sobel(img,cv2.CV_32F,0,1,ksize=3)

#the Gaussian kernel, taken from example program 9
halfSize = 3 * sigma
maskSize = 2 * halfSize + 1
gKern = np.ones((maskSize,maskSize)) / (float)( 2 * np.pi * (sigma**2))
xyRange = np.arange(-halfSize, halfSize+1)
xx, yy = np.meshgrid(xyRange, xyRange)
x2y2 = (xx**2 + yy**2)
exp_part = np.exp(-(x2y2/(2.0*(sigma**2))))
gKern = gKern * exp_part

#convolve the two kernels
Gx = cv2.filter2D(sobelx, -1, gKern)

IxN = my_Normalize(Gx)
cv2.imshow('Gx Normalized', IxN)
return Gx, Ix

img1 = cv2.imread('TestImg1.jpg')
Ix, Iy = my_DerivativesOfGaussian(img1, 1)

cv2.waitKey(0)
cv2.destroyAllWindows()

Here is a side-by-side comparison of the result I was given for a test image (left image), compared with the result I obtained from the above execution (right image):

Where I am confused:

The resulting image that was obtained is close to the target image, except blurred slightly. The blurring leads me to believe there's an error with the Gaussian kernel. I also obtain a drastically different image when I switch the order of the Gaussian and Sobel kernels in the cv2.filter2D() function. However, this is counter-intuitive for me, as the association property states that convolution should be equivalent, regardless of the order in which the convolution occurs.

Specific Questions:

Is the cv2.filter2D() function able to convolve two kernels?
Assuming the answer to ^ is false, do you know of a function that is able to perform such an operation?
Where does my logic/understanding begin to fail?  (i.e. what topics/concepts am I completely not understanding?)

Thanks so much for your time, any help is appreciated.  This is the first time I have been absolutely stumped on an assignment, and can not think of a way to solve this problem without brute-forcing the convolution formula via nested loops.

Here are the images supplied for verifying correctness of the convolution, in their native formats:

Base input test image

Normalized Ix

Normalized Iy

To answer your specifiy questions:

  Is the cv2.filter2D() function able to convolve two kernels?

Yes, but by default, it actually computes the correlation, not the convolution. However, the only difference is that the filter kernel (your second kernel) needs to be flipped. See OpenCV documentation for filter2D.

  Where does my logic/understanding begin to fail? (i.e. what topics/concepts am I completely not understanding?)

I think one more issue is the size of the resulting kernel.
The filter2D function returns an "output image of the same size and the same number of channels as src". However, the combined kernel increases is size (the combined size should be f1-size + f2-size - 1). So before you apply filter2D, you should add a black border to the first kernel.


What happens if we change the limits of integral in Fourier transform? By definition of Fourier transform

$$X(\omega)=\int_{-\infty}^\infty x(t) e^{-j\omega t} dt $$

Now what will happen to the answer of transform  for example in case of $x(t)= \cos(\omega_0 t)$  if limit is $0$ to $A$ instead of $-\infty$ to $\infty$?

For $x(t)=\cos(\omega_0 t)$ its fourier transform is given by $ X(\omega)= \pi[\delta(\omega-\omega_0) + \delta(\omega+\omega_o)]$

so if the limit is changed will it effect the answer?

Yes, it will affect the answer. What you're suggesting is known as the short-time Fourier transform. In the sinusoidal case that you proposed, you will observe spectral leakage, as the truncation of the integral limits is equivalent to multiplication of the sinusoid by a rectangular window function. This multiplication in the time domain maps to convolution in the frequency domain. The Fourier transform of a rectangular window is a sinc function, so the convolution will yield two sinc functions centered at the locations of the impulses in your original answer.


Can FFT convolution be faster than direct convolution for signals of large sizes? Let's say I have a 1D signal of size $N$ and am trying to filter it with a 10-tap FIR filter mask. When $N$ is large, the number of multiply-accumulates would approximately equal
$$2 \times 10 \times N = 20N\quad\text{floating point operations}$$
This convolution algorithm is thus $\mathcal O(N)$ since the computation scales linearly with the input data size.
However, when would it make sense to use an FFT convolution to find this convolution? If an FFT is $\mathcal O(N\log N)$, wouldn't this mean FFT convolution does not scale as well as direct convolution? If this is the case, how can FFT convolution be faster than direct convolution?

There is probably a bit of a misconception here. In many application the signal is  runnning all the time or is VERY long: modem, audio stream, video etc. In this case you can't really define the "length" of the signal. The relevant metric is here "number of operations per input sample" not the "total  number of operations". If you watch a streaming movie, it's obvious that a two hour movie takes twice the operations of a one hour movie. But that's okay since you have twice the time to process it. The key metric is if your processor can keep up with the stream , so that's why "ops per sample" or "ops per second" are used.
In this case, you don't do an FFT over the entire input signal but you apply "overlap add", i.e. you chop the input signals into blocks that are about the same size as the filter length and then you process one block at a time.
So the effectiveness of the algorithm is only a function of the filter length. And it increases linearly for direct convoltuion but only logarithmically for Overlap Add (which is FFT based). The "break even" point is typically at "many tens of samples".
We can run some numbers assuming that the filter length is $N$ and lets also assume its a power of 2.
An FFT of length $M$ has  $log_2(M) \cdot M/2$ butterflies, each with one complex multiplies and two complex adds. That's a total of 10 real operations. So the price for an FFT is roughly $5 \cdot log_2(M) \cdot M/2$
Overlap add requires zero padding to 2N and that's our FFT size. Hence overlap add costs

Forward FFT: $5 \cdot log_2(2N) \cdot 2N/2$
Spectral Multiply: $3*N$  (N/2 complex multiplies)
Inverse FFT: $5 \cdot log_2(2N) \cdot 2N/2$
Overlap adding: $N$
Total: $N\cdot (10 \cdot log_2(N) + 10 + 3 + 1)$

So overlap add cost per sample is roughly $ 10 \cdot log_2(N) + 14$ as compared to $2N$ for direct convolution. The break even point is a filter length of 32.
DISCLAIMER: I made a lot of simplying assumptions and in a practical implementation the break even point can vary wildly based on processor architecture and smartness of the algorithm. The main cost is in the FFT and for a single channel real system you would deploy a core complex FFT of N/2. I also skipped over bit reverse ordering and other overheads.


What is the best way to detect rectangular pulse for AFSK demodulation? I've got a WAV file with a rectangular pulses coding zeros and ones. By specification, a 1 bit is a 320 microseconds wide rectangle, and a 0 bit is a 640 microseconds.

So basically the idea is to detect "zero points" where the signal comes over the time axis. When I get it, further processing is trivial.

The naive approach is to iterate samples in the signal and find points where sample * sampleNext < 0.

For low-noise files with this approach produces rectangles quite close to the specified ones, but for high-noise ones I get all sorts of rectangles from 20 to 1400 microseconds in width.

Is there a better way to detect zero points in a noisy rectangular pulse?

How about you take a short time Fourier transform with window length equal to the granularity with which you want to locate the zero crossing. At a zero crossing the signal will have high frequency components. You could detect those by running a short time Fourier transform on the incoming samples to identify the location of zeros. For ex: if you window 5ms of data and take FFT and see if there are high frequency components of significant magnitude then this would mean that this window contains the zero crossing. If you want, you can go even narrower like 1ms, but you be would loosing on frequency resolution. So, pick an optimum resolution in frqeuency and required zero crossing accuracy in time and window the incoming signal and keep taking the FFTs. You could go for rectangular windows.

Even if noise is high it will be constant in spectrum (atleast close to it as is theoretically for white nosie) and this should not be a big trouble in FFT domain, unless it is higher compared to signal power itslef.

For just 13 sample long pulse durations use the following emperical algorithm

Maintain two numerical derivate with a reference point x samples apart. Fir ex: consider "sample number" $N_o$, take two derivates between $N_o$, $N_o -x$ call this gradient 1, and second derivative between $N_o$, $N_o + x$, call this gradient 2, take the magnitude of these two gradients, if they are greater than a threshold, usually gradient is close to 90 degree near the  zero crossing, then the point $N_o$ is a good approximation of zero crossing.  $x$, and the threshold can be tuned based on pulse length.


People counting implementation using OpenCV I need to implement people counting(head count) into one of the product.
I am able to detect faces and count them but only when the persons are actually looking into the camera, if they are not looking into it they will not be detected.

So, what I want is to an idea about detecting human rather than detecting only face.

I have gone through some algorithm like voila jones etc..but did not get anything useful.

It would be great if someone can share their experience or any code snippet to do so or any algorithm to help me out here.

This is an ongoing research problem, which is much more difficult to solve than face detection due to a gigantic variety of poses, angles, and other factors that play a role. This is usually done with image feature extraction on a large set of images followed by a supervised learning algorithm. In short, there's no unique answer, and there's definitely not a simple "code snippet" one could use for this task.

If you're implementing a specific algorithm, you could probably get better help here, but not at this point.


What is different between gradient and imgradient functions in MATLAB I would like to compute the gradient of an image I. I have the following two options :

[Gx, Gy] = gradient(I);
g = sqrt(Gx.^2+Gy.^2);

or

[g,~] = imgradient(I, 'sobel');

My questions are

What is the gradient method used in the first option ?
What are the benefits of using Sobel method to estimate the gradient ?

Question 1:
gradient returns both estimated horizontal $\left(\frac{\partial I}{\partial x}\right)$ and vertical components $\left(\frac{\partial I}{\partial y}\right)$ of the gradient vector, which in your example code are Gx and Gy. To do so, gradient computes for each point $X(i)$ the half difference of $X(i-1)$ and $X(i+1)$, which is the most basic gradient estimation method (it's called the central difference method).
Then, computing

g = sqrt(Gx.^2+Gy.^2);

gives you the gradient magnitude at each point.

imgradient differs since it directly outputs the gradient's magnitude (g) and direction (which you suppressed using ~) at each point. In addition, it computes the gradient using a two-dimensional method applying a sobel kernel to each point and its neighbors. Finally, this function allows you to specify which kernel to use to compute the numerical gradient if you don't want to use the Sobel estimator. You can for example choose to estimate the gradient with the central difference method as in gradient.

Question 2 :
Sobel method is just the simplest contour detection method providing good results for a low computation cost. But there are plenty of others, that are used depending on the application I guess.
Compared to the central difference method, it achieves better noise rejection using a different kernel with coefficients designed to do so. See as an example the kernel used for horizontal gradient estimation :
$$
\begin{matrix}
-1 & 0 & 1\\
-2 & 0 & 2\\
-1 & 0 & 1
\end{matrix}
$$
The considered point's row is emphasized (factor 2) compared to it's neighbors (factor 1). While the central difference method does not make that difference and weights all points with the same coefficient :
$$
\begin{matrix}
-1 & 0 & 1\\
-1 & 0 & 1\\
-1 & 0 & 1
\end{matrix}
$$


Explaining / Communicating spectra: Generating Visualizations I'm very often in the situation that I need to explain / visualize spectrum of signals or frequency behaviour of systems (esp. filters). My favourite, very "textbook" styled way of communicating these would be, in simpler cases:

¹

Typically, I'd have to visualize something like aliasing, effects of suppressing carriers etc. This also leads to less simple situations, like

¹

Now, drawing these kinds of diagrams is really a burden if done with my vector program of choice (that being Inkscape) – I want to parametrically specify the spectrum, not manually shift, copy, mirror and repeat lines!

So:

How do you approach theses kind of visualizations? Vector graphics programs & patience? Pen & paper? Or do you have a private set of scripts that e.g. draw the diagram with the help of matplotlib or tikz?

¹ Bellanger: Digital Processing of Signals: Theory and Practice. 2nd Edition, Wiley 1988. P. 17 & p. 18.

Personally, I value aesthetics and consistency in my documents, so I use Latex whenever I can. I use Inkscape if I'm in a hurry, though.

I have found that the pgfplots package (built on top of tikz) reduces the time needed to code a plot substantially, especially once you get the hang of it. As a simple example, the following filter response:

was generated with the following Latex code:

\documentclass{standalone}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis} [
        height=5cm, width=7cm, xlabel={$f$}, ylabel={Filter A},
        xmin=0, xmax=300, ymin=0, ymax=1.6, no markers,
        xtick={100,175,250}, ytick={0.25,0.5,0.75,1},
        axis lines=middle
    ]
    \addplot coordinates {(0,1) (100,1) (250,0)};
    \addplot[black,dotted] coordinates {(0,0.5) (175,0.5) (175,0)};
\end{axis}
\end{tikzpicture}
\end{document}

Excluding boilerplate (which can be reused between plots), this is just two lines of "real" code (the \addplot lines). The benefits are the usual latex quality, plus font/style/size consistency with the rest of the document.

As a further example I tried a (simplified) version of the last plot in the question:

which requires just a bit more code:

\documentclass{standalone}
\usepackage{pgfplots}
\usepgfplotslibrary{fillbetween}
\begin{document}
\begin{tikzpicture}
\begin{axis} [
        height=5cm,width=9cm,xlabel={$f$},axis lines=middle,
        xmin=-120,xmax=200,ymin=0,ymax=1.8,no markers,
        xtick={80,160},xticklabels={$\frac{f_s}{2}$,$f_s=\frac{1}{T}$},
        ytick={0},yticklabels={},ylabel={$S_s(f)$}
 ]
    \addplot[black] coordinates {(-100,0) (-20,1) (0,0)};
    \addplot[black] coordinates {(0,0) (20,1) (100,0)};
    \addplot[black,name path=A] coordinates {(60,0) (140,1) (160,0)};
    \addplot[black,dashed,name path=B] coordinates {(80,0) (80,1.5)};
    \addplot[blue] fill between[of=A and B,soft clip={domain=60:80}];
\end{axis}
\end{tikzpicture}
\end{document}

It is true that generating these plots takes time. When working on anything but the most trivial plot, I constantly need to refer to the manuals and to tex.stackexchange.com. I think the examples above show that the pgfplots package removes some of the pain; once you get practice with the most common commands, it shouldn't take too long to create nice plots.

A final point: while I don't (yet) have a program to automatically generate spectrum plots, it shouldn't be too hard to do. I do have a script to generate pgfplots scripts from Matlab/Octave, which you can see here. This script takes two or more Matlab arrays and a plot description, writes a latex file, and runs pdflatex on it. While laborious, coding it wasn't particularly difficult. It should be possible to modify it to take spectrum descriptions (maybe in the edges, amplitudes format of the filter design commands) and generate the plots automatically.


Using overlap add for drawing spectrogram I'm confused with overlap-add (OLA) in the following two scenarios

filtering/convolution an input signal;
applying an STFT on input signal only and draw its spectrogram.

When applying OLA to filtering/convolution, I'd take overlapped signal blocks, window them, take FFT of each windowed block, multiply filter core in frequency domain, IFFT the block, then apply  OLA to the current block with buffered previous block in time domain.

When drawing spectrogram of the signal, I assume I must take overlapped signal blocks, window them, take FFT of each windowed block. Now I'm not sure how I would draw the spectrogram with the magnitude response with the block.

Which way of the following is correct?

Apply window function to each overlapped window, and then simply draw the FFT magnitude response of each window, making the windows next to each other as if the spectrum of each window corresponds to the center time or left edge of the window.
Apply window function to each overlapped window, add up the windowed time-domain signal between adjacent blocks. Then take FFT of the current block and draw it. This requires caching one signal buffer.

The first approach seems to give the wrong timeline, if we consider the block signal starts at its left edge and ends at its right edge, because it's overlapped with the next block, unless we say the resulting frequency response represents the non-overlapped portion of the window only.

The second approach seems to counteract the windowing and thus beats the purpose, because adding the overlapped portion in time domain will cancel out windowing on the overlap.

I'm confused which way would be the correct OLA for deriving the spectrogram of the running signal.

Apply window function to each overlapped window, and then simply draw
  the FFT magnitude response of each window as if the spectrum
  corresponds to the center time of the window.

Your approach is partially correct. But in the second half of sentence is not how you draw a spectrogram. You shouldn't align the FFT output in the same direction as time. It should be perpendicular direction as time. After you obtain FFT Magnitude, look at it from "above". Means the colour grades that you give depends on the magnitude value. This will be line corresponding to first window. Similarly, for second window, you place it next to it (since you took it after certain interval of time). An example below shows simple spectrogram having 3-bin FFT with middle bin having higher value compared to other bins.


How does PSOLA pitch-shift work in real-time? How does the PSOLA pitch-shifting algorithm manage to remove or add audio frames in real-time without creating gaps?

Wouldn't extending the "windows" (it's like in FFT windowing) to longer than the audio buffer create asynchronity in the buffering? Similarly compressing the windows would lead to the need to "read new audio in advance" in order to not leave silence.

Or what's the idea?

I guess there's some clever math (ratios?) in order to know how to fill and remove windows so that they still match up reasonably in the overlap-add process.

In PSOLA, pitch estimation is used to space out or overlap segments by integer multiples of the estimated pitch period in time.  Previous frames (or pitch period length sub-segments thereof) can be duplicated (as needed) in order to not leave any gaps between frames.  Segments can then be concatenated or cross-faded.  This can lengthen or shorted the composited result, and thus change the time duration.  Then this intermediate result can resampled to change the pitch (and possibly undo some or all of the time duration change), and to provide the proper number of samples needed for the output sample rate.

There will be some jitter, as frames are usually only shifted by a minimum of one whole pitch period, not by fractions thereof.  Buffering may be required to cover this time jitter, the amount related to the longest period (lowest pitch) expected to be handled.

The "clever math" might be in the choice and implementation of an accurate enough but low latency pitch estimation algorithm.  If the pitch cannot be estimated (due to consonants or polyphony, etc.), apparent quality may suffer.


Is this system time invariant? I've been working at this problem for a while now, and can't seem to come to a solid conclusion - is this system time invariant?

$y(t) = \int_{-\infty}^{t} e^{-9(t-\tau)} x(\tau)d\tau $

My reasoning for the system being time invariant is that it is clearly a convolution integral. Time shifting the input will shift the output.

But, it is not a normal convolution integral. The top bound is only to $t$ rather than to $\infty$. This is all well and good if the system is causal, because for any $\tau$ greater than $t$, the output would be zero. But because there is no unit step $u(t)$ clearly defined in the transfer function (like $e^{-9(t-\tau)}u^\tau(t)$), one cannot make a conclusion on the causality of the system. So the question at hand is, does shifting the input also shift the bounds of the integral?

So what do you think? Can anyone prove if the system is time invariant as given?

Hint

You are on the right track. You note that the integral is

  not a normal convolution integral

and

  there is no unit step $u(t)$ clearly defined in the transfer function.

So, can you introduce a unit step yourself and change the limits of integration to make it look like a normal convolution integral?

STOP! Full answer below.

Let $\mathbf{1}_{A}$ be the indicator function which is 1 on the set $A$ and 0
everywhere else.
\begin{eqnarray}
y(t) &=& \int_{-\infty}^t e^{-9(t-\tau)} x(\tau) d\tau \\
&=& \int_{-\infty}^{\infty} e^{-9(t-\tau)} x(\tau) \mathbf{1}_{\{\tau \leq t\}} d\tau \\
&=& \int_{-\infty}^{\infty} e^{-9(t-\tau)} u(t-\tau) x(\tau) d\tau
\end{eqnarray}

This is an LTI system with transfer function $e^{-9t}u(t)$.


Soft question: Why do we need to reconstruct a signal? I am currently reviewing the Nyquist criterion which says I need to sample the signal at 2 times the maximum frequency in order to avoid aliasing in the reconstructed signal.

My question is that in all the applications I've done, I rarely ever remember having to reconstruct the signal after sampling it.

For example, in controlling a digital controller, all I have to do is to sample the signal and send the sampled signal to the controller. End of story, no reconstruction required.

In other applications, I sample the signal to analyze its frequency content. After looking at the frequency distribution, the signal is completely discarded. No reconstruction required.

Can someone make it clear why do we need to respect the Nyquist criterion if no reconstruction is needed (so no worry about aliasing effect)?

The concept of reconstruction has nothing to do with the application, rather it has to do with the question: did I get the same signal that is really there. If you cannot recreate the signal back, that means the conversion process is loosing/modifying underlying information, which in most cases you do not want to happen.

So the confidence on the Fourier Transform, or to use the digitized data in the controller, comes from the fact that the loss of information is within the tolerance limits of your application. You can measure the quality of your digitized signal, by actually re-synthesizing the signal and comparing it to the original signal.


Reversing the order of up and downsampling When we are dealing with the problem of changing the sampling rate by a noninteger factor like L/M, we first need to the do the upsampling, and then the downsampling.

My question is that under which condition can we reverse this order?

The problem with downsampling is that it can be lossy -- since you're reducing the sampling rate, you can introduce aliasing. So, you can reverse the order whenever downsampling does not result in aliasing.

For example, say your discrete-time signal $x[n]$ contains energy in frequencies up to $f_N/3$, where $f_N$ is the Nyquist frequency. Then, downsampling by a factor of 2 will not affect the signal: the "new" Nyquist frequency is $f_N/2$, so there is no aliasing.

See also Which order to perform downsampling and filtering?.


Minimum number of points to specify a sine wave with a given accuracy I am generating a sine wave from a DAC at a constant frequency. If each full cycle is specified by N data points what is the relationship between the accuracy of the rendered sine wave and N?

N=1 will obviously be pulses. N=4 (peaks and zero crossings) will yield a triangular wave of sorts. Is there some equation that shows how good a sine wave approximation is being created? [Obviously I can put it through a post DAC filter to smooth it etc]

The kind of answer I am looking for is: "Specifying a sine wave with N points yields 1:X accuracy"

I think what you are asking is "what is the sample rate required to sample a sine wave for a given precision". This is basically what the sampling theorem is all about (See for example http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem). The surprising answer is (more or less): "more than two". Once you meet the sampling criteria you are fine and
the accuracy does NOT increase by sampling more often.

In your case, 3 points per period will work fine (provided your DAC has a proper anti-aliasing filter).


Difference between Sub Sampling and Down Scaling of Images I know that there are many ways to upscale (interpolate) an image using bilinear, bicubic, sinc... Somehow, these same algorithms can also be used to downscale an image. But when it comes to subsampling, I've come across only two methods: Gaussian-blur then discard some rows and columns, or block-averaging.

Is there such thing as bilinear, bicubic,... subsampling?

Instead of Gaussian-blurring then discarding or block-averaging, can/should I do Gaussian-blurring then block-averaging? or bilinear, bicubic,...subsampling if there is such thing?

What I'm doing is trying to make image Gaussian and Laplacian pyramids; but I'm not sure how to decimate and interpolate the image. (Some Python code with Scipy/Numpy would really help.)

Thank you very much.

Have you looked to the pyUp and pyDown functions in opencv? They upscale and downscale images using a gaussian pyramids.

opencv image filtering page


PortAudio: how to generate multiple sounds I have simple wave generator in PortAudio. I'm creating single sine wave, and then I use a lookup sine to play it. Now I've stucked here. I have no idea how I can add another sine sound. Is it possible to somehow add new channel with another sound?

Here is my code:

soundengine.cpp

#include "soundengine.h"
#include <QDebug>

void SoundEngine::run()
{
    PaStreamParameters outputParameters;
    int i;
    double t;
    turnOFF = false;
    isPressed = false;

    static unsigned long n=0;
    for( i=0; i<TABLE_SIZE; i++, n++ )
    {
        t = (double)i/(double)SAMPLE_RATE;
        data.sine[i] = 0.80f * (float) sin( ((double)i/(double)TABLE_SIZE) * M_PI * 2. );
    }
    data.phase = 0;
    data.phase_increment = CalcPhaseIncrement(0);

    err = Pa_Initialize();
    if(err != paNoError) qDebug()<<"Błąd przy inicjalizacji strumienia:"<<Pa_GetErrorText(err);

    outputParameters.device = Pa_GetDefaultOutputDevice(); /* default output device */
    if (outputParameters.device == paNoDevice) qDebug()<<"Błąd: Brak domyślnego urządzenia wyjścia!";

    outputParameters.channelCount = 2;       /* stereo output */
    outputParameters.sampleFormat = paFloat32; /* 32 bit floating point output */
    outputParameters.suggestedLatency = Pa_GetDeviceInfo( outputParameters.device )->defaultLowOutputLatency;
    outputParameters.hostApiSpecificStreamInfo = NULL;

    err = Pa_OpenStream(
              &stream,
              NULL, /* no input */
              &outputParameters,
              SAMPLE_RATE,
              FRAMES_PER_BUFFER,
              paClipOff|paDitherOff,      /*paNoFlag we won't output out of range samples so don't bother clipping them */
              patestCallback,
              &data );
    if(err != paNoError) qDebug()<<"Błąd przy otwieraniu strumienia:"<<Pa_GetErrorText(err);
    err = Pa_StartStream( stream );
    if(err != paNoError) qDebug()<<"Błąd przy starcie strumienia:"<<Pa_GetErrorText(err);

    while (turnOFF == false) {
        Pa_Sleep(500);
    }

    err = Pa_StopStream( stream );
    if(err != paNoError) qDebug()<<"Błąd przy zatrzymywaniu strumienia:"<<Pa_GetErrorText(err);
    err = Pa_CloseStream( stream );
    if(err != paNoError) qDebug()<<"Błąd przy zamykaniu strumienia:"<<Pa_GetErrorText(err);
    Pa_Terminate();
}

float SoundEngine::LookupSine( paTestData *data, float phase )
{
    float fIndex = phase*TABLE_SIZE;
    int   index = (int) fIndex;
    float fract = fIndex - index;
    float lo = data->sine[index];
    float hi = data->sine[index+1];
    float val = lo + fract*(hi-lo);
    return val;
}

int SoundEngine::patestCallback(const void *inputBuffer, void *outputBuffer, unsigned long framesPerBuffer, const PaStreamCallbackTimeInfo *timeInfo, PaStreamCallbackFlags statusFlags, void *userData)
{
    paTestData *callData = (paTestData*)userData;
    float *out = (float*)outputBuffer;
    float sample;
    unsigned long i;

    (void) timeInfo; /* Prevent unused variable warnings. */
    (void) statusFlags;
    (void) inputBuffer;

    for( i=0; i<framesPerBuffer; i++ )
    {
            sample = LookupSine(callData, callData->phase);
            *out++ = sample;  /* left */
            *out++ = sample;  /* right */

        callData->phase += callData->phase_increment;
        if( callData->phase >= 1.0f ) callData->phase -= 1.0f;
    }

    return paContinue;
}

void SoundEngine::generateSine()
{
    if(isPressed == false)
    {
        data.phase_increment = CalcPhaseIncrement(220.0);
        isPressed = true;
    }

}

void SoundEngine::removeSine()
{
            data.phase_increment = CalcPhaseIncrement(0.0);
            isPressed = false;
}

To play sound I use generateSine() method. Maybe is there better way to play sound at given frequency, using one method(like generateSine() in my example). I would like to play two sounds by using generateSine() method twice.

The beauty of PCM (pulse-coded modulation, i.e. time samples of audio) is that the mathematical sum of two signals sounds like the acoustical superposition.

So, just add your two sines sample-wise.


When to use the DTFT vs the DFT (and their inverses) in analysis? In many of my readings, whenever some author speaks about working in the frequency (transform) domain (of a digital signal), they often times take the DFT, or the DTFT, (and of course their corresponding inverses). Different authors will tend to work with one or the other.

I have not been able to really ascertain a particular pattern regarding this. In that, why would you pick the DTFT over the DFT or vice versa in explaining algorithms? Where does one help you over the other?

DFT and DTFT are obviously similar as they both generate the fourier spectrum of time-discrete signals. However, while the DTFT is defined to process an infinitely long signal (sum from -infinity to infinity), the DFT is defined to process a periodic signal (the periodic part being of finite length).

We know that the number of frequency bins in your spectrum is always equal to the number of samples processed, so this also gives a difference in the spectrums they produce: the DFT spectrum is discrete while the DTFT spectrum is continuous (but both are periodic with respect to the Nyquist frequency).

Since it is impossible to process an infinite number of samples the DTFT is of less importance for actual computational processing; it mainly exists for analytical purposes.

The DFT however, with its finite input vector length, is perfectly suitable for processing. The fact that the input signal is supposed to be an excerpt of a periodic signal however is disregarded most of the time: When you transform a DFT-spectrum back to the time-domain you will get the same signal of wich you calculated the spectrum in the first place.

So while it does not matter for the computations you should note that what you are seeing there is not the actual spectrum of your signal. It is the spectrum of a theoretical signal that you would get if you periodically repeated the input vector.

So I would assume in the literature you were mentioning, every time it is important that the spectrum you're working with is actually the spectrum and disregarding the computation side of things, the author would pick the DTFT.


How to smoothly adjust sine amplitude I move a device by sending a sine wave proportional to the displacement. I would like to modify the amplitude and mean value dynamically.

For example at the beginning I would like to go from :

Signal(t) = 1 * sin(t) + 1

to:

Signal(t) = 0.95 *sin(t) + 0.95

My issue is that I don't know how to do that without having a big step in my signal, which is unacceptable for my application.
For now, all I have is something that change the coefficient on a "long" time, by small step (1>0.99>0.98>0.97>0.96>0.95) so that the "steps" in the output signal remain small, but it is not very clean.

There's no problem letting the coefficients be a slowly varying function of time,

$$S(t)=a(t)\sin(t)+b(t),$$

for example a linear ramp
$$a(t)=(a_1-a_0)\frac{t-t_0}{t_1-t_0}+a_0$$ for $t\in[t_0,t_1]$ and similar for $b(t)$.


Simulation of a Frequency ramp This might be an easy one!

Please preface your answer with spoiler notation by typing the following two characters first ">!"

Given the following equation that describes a cosine wave as a function of time:

$$m(t) = \cos(2\pi f t) $$

We know that $t$ is the independent variable for time in seconds, and $f$ as a constant value will be the frequency of the cosine in $\textrm{Hz}$.

What if f is a linear ramp over time?

Given a time domain waveform $m$ with a linearly increasing frequency over time, that starts at $t=0$ with frequency $F_1\textrm{ Hz}$ and ends with frequency $F_2\textrm{ Hz}$:

Provide the expression for $f(t)$ in terms of $F_1$ and $F_2$ and $T$, that would create a function $m(t)$ whose frequency starts at $F_1$ and ends at $F_2$, with the frequency linearly increasing over some arbitrary time period $T$.

The actual spectrum, as in any spreading effect due to the ramp rate, is of no interest to this solution.

Hint:

   The answer is a simple linear equation ($y = mx+b$), with a $y$-intercept $b$ that is a constant dependent only on $F_1$ and value $y$ at time $T$ is another constant dependent only on $F_2$  --- what are those two constants with relation to $F_1$ and $F_2$?

I waited a bit to see if someone else takes the challenge, but since there are no answers yet, I'm providing mine now.

   For a given signal $m(t)=\cos(\phi(t))$ the instantaneous frequency is given by $$\nu(t)=\frac{\phi'(t)}{2\pi}\tag{1}\\$$One pitfall here is that if we write $m(t)=\cos[2\pi f(t)t]$, $f(t)$ is generally not the instantaneous frequency (unless $f(t)$ is constant). We want the following instantaneous frequency:$\nu(t)=\begin{cases}F_1,&t<t_1\\F_1+\frac{F_2-F_1}{t_2-t_1}(t-t_1),&t_1<t<t_2\\F_2,&t>t_2\end{cases}\tag{2}\\$From $(1)$ and $(2)$, in the time interval $[t_1,t_2]$ the phase $\phi(t)$ must be chosen as $$\frac{\phi(t)}{2\pi}=F_1t+\frac{F_2-F_1}{t_2-t_1}\left(\frac{t^2}{2}-t_1t\right)+C\tag{3}$$Choosing $C=0$, and since $\phi(t)/(2\pi)=f(t)\cdot t$, we get from $(3)$ $$\begin{align}f(t)&=F_1+\frac{F_2-F_1}{t_2-t_1}\left(\frac{t}{2}-t_1\right)\\&=F_1-\frac{\Delta F}{\Delta t}t_1+\frac12\frac{\Delta F}{\Delta t}t,\quad t_1<t<t_2\tag{4}\end{align}\\$$ with $\Delta F=F_2-F_1$ and $\Delta t=t_2-t_1$.


Continuous vs discrete signal energy I am looking to calculate the signal energy of real, sampled acoustic data. According to this source, the energy of a continuous signal is:
$$\tag{1}\hspace{1cm} E_x = \int_{-\infty}^{+\infty}|x(t)|^2dt $$
This is intuitive as we are often taught that the energy is given by the area under the curve, for a power curve. Certainly the equation 1 can easily be approximated by the Riemann Sum which still retains the differential
$$\tag{3}\hspace{1cm} S = \sum_{i=1}^{n} f(x_i)\Delta x_i $$

However, according to this source, the discrete energy calculation is given below as:
$$\tag{2}\hspace{1cm} E_x= \sum_{-\infty}^{+\infty} |x(n)|^2$$

My question is why do we drop the differential ($dt$)?

I think you are correct. People are being fast and loose with the expression in your Eq. (2), but it captures the behaviour of the energy of the signal up to a constant $T_s$ factor (the sampling period), which is maybe why they do it.

As you say, the signal energy is given by
$$
E_s = \int_{-\infty}^{+\infty} |x(t)|^2\;dt \tag{1}
$$
If you want to approximate this using a Riemann sum for the discrete case, the equivalent expression (depending on your choice of indexing) is
$$
E_s \approx \sum_{n=1}^{\infty} x[n]\; x^*[n]\cdot T_s \tag{2}
$$

You can see immediately that there must be a factor of $T_s$ in there otherwise the dimensions of the units aren't consistent. The correct units of signal energy as given by Eq. (1) are $[\textrm{signal}^2\cdot\textrm{time}]$. If you drop the $T_s$ then the units of your Eq. (2) would just be $[\textrm{signal}^2]$, which are NOT the correct units for signal energy. See this related question for the units.

EDIT---------------------------------

The OP cited the wiki reference as the source of the confusion, specifically with the discrete version of Parseval's theorem. Parseval's theorem in continuous time says that the energy in the time domain and frequency domain are equal, and is written as

$$
\int_{-\infty}^{+\infty} |x(t)|^2\;dt = \int_{-\infty}^{+\infty} |X(f)|^2\;df \tag{3}
$$
where $x(t)$ is the signal, and $X(f)$ is its Fourier transform, defined by
$$
X(f) = \int_{-\infty}^{+\infty} x(t) e^{-2\pi i f t}\; dt. \tag{4}
$$
You can verify together with Eq. (4) that the units of both sides of Eq.(3) are $[\textrm{signal}^2\cdot\textrm{time}]$, as they should be for signal energy. So far so good.

You can get a discrete version of Parseval's theorem by discretising both sides to get
$$
\sum_{n=0}^{N-1} |x[n]|^2\cdot T_s = \sum_{k=0}^{N-1} |X[k]|^2\cdot \Delta f \tag{5}
$$
where $T_s$ is the sampling period, and $\Delta f$ is the separation of frequency bins (I have assumed that the length of the vector $X[k]$ and $x[n]$ are both equal to $N$ - this would be different if you had zero-padding, but remember that the extra zeroes would not contribute any extra signal energy). The Fourier transform $X[k]$ is the discrete version of $X(f)$, and is given by
$$
\begin{align}
X[k] &= \sum_{n=0}^{N-1} x[n] \exp \Big(  -2\pi i \frac{kn}{N} \Big) \cdot T_s \tag{6}\\
&= X_{\textrm{DFT}}[k]\cdot T_s \tag{7}
\end{align}
$$
and we recognise $X_{\textrm{DFT}}[k]$ as being the result of the discrete Fourier transform (DFT) operation - this you might obtain directly from the output of the FFT function, for example in MATLAB. There is an extra factor of $T_s$ in the expression for $X[k]$, because it is an approximation of the continuous Fourier transform from Eq. (4) - again needed for the units to add up.

Now we can insert Eq. (7) into Eq. (5), and use the fact that $\Delta f = F_s/N = 1/(NT_s)$, where $F_s$ is the sampling frequency, to obtain
$$
\begin{align}
\sum_{n=0}^{N-1} |x[n]|^2\cdot T_s &= \sum_{k=0}^{N-1} |X_{\textrm{DFT}}[k]\cdot T_s|^2\cdot \Delta f \tag{8}\\
&= \sum_{k=0}^{N-1} |X_{\textrm{DFT}}[k]|^2\cdot \frac{T_s^2}{NT_s}\tag{9}\\
\implies \sum_{n=0}^{N-1} |x[n]|^2 &= \frac{1}{N}\sum_{k=0}^{N-1} |X_{\textrm{DFT}}[k]|^2\;,\tag{10}
\end{align}
$$
where the factors of $T_s$ have cancelled. This is the expression for Parseval's theorem given in the wikipedia link you cited, and is still valid. However, since we killed the factors of $T_s$, the units are now not signal energy units, and I personally wouldn't call the left hand side of Eq. (10) signal energy anymore (actually, it has units of signal power instead).

Hope that helps clear up a bit of the confusion. I always find it helpful to check the units at each stage as you go along, to keep track of things.


Is it possible to combine two or more Biquads? I've wrote a Equalizer, and now it apply the Biquads in serie (cascade), I was just thinking if i can combine the biquads it will more efficient than applying one-per-one.

Is that possible?

Possible, yes.  But higher-order IIR filters tend to be less stable and more susceptible to various numerical problems (lost precision, underflow, overflow, etc.), compared to a cascade or bank of suitable partitioned biquad filters.


Finding maxima and minima of a continuous time signal by using frequency method? I have a PPG signal and I wanted to write a code that gives the maxima and minima value of signal by using frequency domain method. Any suggestions?

Global maximum

Representation of a signal from DFT components is:

$$f[n] = \frac{1}{N} \sum_{k=0}^{N-1} F[k] e^{i 2\pi k n /N}$$

Let $\theta = 2\pi n /N$

$$f[N\theta/(2\pi)] = \frac{1}{N} \sum_{k=0}^{N-1} F[k] e^{i k \theta}$$

This is a trigonometric polynomial with coefficients $F[k]/N$. Finding the maximum is the same as for a normal polynomial, at the maximum the derivative is 0. Thus you calculate the derivative polynomial,

$$\frac{1}{N} \sum_{k=0}^{N-1} i k F[k] e^{i k \theta}$$

Use root finding to obtain $N-1$ possible values of $\theta$ (the roots), substitute them back into the original polynomial and choose the one that corresponds to the maximum. Then $n$ can be found as $n = N\theta/(2\pi)$.

Global minimum

Same as the above except you choose the value of the original polynomial that is a minimum

Local maxima and minima

Find the $N-1$ roots of the polynomial as before. Discard any points which do not correspond to derivative of zero. The remaining points are local maxima or minima. At the local maxima the second derivative is negative, at local minima it is positive.


What are the advantages and disadvantages separability of Gaussian 2D filter? I know separability of a Gaussian filter enhances the computational complexity from $\mathcal{O}(L^2*N*M) \to \mathcal{O}((L*N*M)$. How this really reflect on my program? I suppose the number of for loops that I need for the convolution less compare to nonseparable case? What real disadvantages ?

Assuming the simplest case with a square image $x[n,m]$ of size $N \times N$ and a square filter kernel $h[n,m]$ of size $M \times M$, the raw 2D convolution to produce the, cropped, output image $y[n,m]$ of size $N \times N$ requires about $N^2 \times M^2$ MACs.

The raw 2D convolution between $x[n,m]$ and the filter $h[n,m]$ is implemented by using two for loops to range through each output pixel and two additional for loops to perform the 2D convolution at that pixel location. Hence a total of 4 nested for-loops are required; a complexity of $\mathcal{O}(N^2M^2)$

When the filter kernel $h[n,m]$ is separable, such as $h[n,m] = f[n]g[m]$, then in this case the convolution between the image $x[n,m]$ and the filter $h[n,m]$ can be performed without a raw 2D convolution sum, by the following approach:

First, perform a 1D convolution between columns of $x[n,:]$ and the 1D filter $f[n]$, which requires about $N \times M$ MACs to complete. This operation should be performed for each column of $x[n,m]$ by proceeding along its horizontal, $N$ many, columns. Hence a total of $N \times M \times N$ MACs will be required to complete the first step to produce the intermediate image $w[n,m]$
Then, apply the similar algorithm, with the filter $g[m]$ and rows of the intermediate image $w[:,m]$, which will require similar number of MACs as $N \times M \times N$.

Hence as a total you will need about $2N^2M$ MACs for the separable implementation of the 2D convolution. The actual number depends on the cropping type you apply. Thus a complexity of about $\mathcal{O}(N^2M)$ is attained.

It can be seen that, to implement the separable convolution algorithm, you will only need 3 nested for loops, which is where the gain comes from.

The main advantage of separable filtering is quite clear; much reduced computational cost. In fact even the 2D-FFT algorithm makes use of it as the 2D-DFT kernel is separable.

Another advantage is that since there are less number of MACs involved to produce the same output compared to the raw 2D convolution, it's less prone to numerical issues and likely to produce more accurate results.

A disadvantage of separation is that it requires an extra ram memory to store the intermediate image, which can be a concern in certain applications.


Demodulation of FSK signal What kind of FSK signal is this and what demodulation technique can I use to demodulate it. Why is the phase of the 1's changing like that and how does that effect how it needs to be demodulated. I will be programming the demodulation in python.

I am trying to low pass at 900Hz and then generate the envelope but I am not getting the desired result.

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
from scipy.fftpack import fft, rfft, rfftfreq, irfft
import scipy.signal.signaltools as sigtool

from scipy.io import wavfile
fs, data = wavfile.read('jx3p.wav')

h = signal.firwin(numtaps = 300, cutoff = 900, fs=fs)
data = signal.lfilter(h, 1, data)

data = np.abs(sigtool.hilbert(data))

import matplotlib.pyplot as plt
plt.plot(data)
plt.show()

UPDATE

Here is the code to multiple the signal by the same signal with a 12 sample delay and then add the resulting signal with a 3 sample delay of itself.

import matplotlib.pyplot as plt
from scipy.fftpack import rfft, rfftfreq
from scipy.io import wavfile
from scipy.signal import blackman

fs, data = wavfile.read('jx3p.wav')

x = 286600
y = x+22050
data = data[x:y]

m = max(abs(data))
data = [d/m for d in data]

delay = 12
data1 = tuple(a*b for a,b in zip(data[delay:],data))

filter_delay = int(delay/4)
data2 = tuple(a+b for a,b in zip(data1[filter_delay:],data1))

plt.plot(data[:300],"r")
plt.plot(data1[:300],"b")
plt.plot(data2[:300],"g")

plt.show()

For this specific waveform as described, the following would demodulate the signal from the Frequency Shift Modulated Input into a square wave output:

This works given the duration of the "0" symbol is 12 samples by using a "delay and multiply" frequency discriminator. The multiply will have a strong double frequency component that needs to be filtered out (and can be used for timing recovery if the timing was not synchronous to the symbol rate as further detailed below), and this is done with a simple delay of 3 samples and add, which would provide a null at twice the frequency of the zero symbol.  The "1" symbol as shown by the OP is a low frequency symbol in the FSK modulation, and observe how in this case it is approximately formed using the same half cycle of the higher frequency "0" symbol followed by a zero magnitude for a full one and a half cycle duration before continuing with the negative half cycle followed by a zero magnitude for a duration for one and a half cycles (thus it is 4 times the duration of the "1" symbol.) This is convenient as the output for the case of demodulating the "0" symbol would be naturally nulled given the block diagram above without the need for subsequent filtering.

Notice with this approach that at the start of each symbol, if the symbol is a 1 the initial first half cycle result will be a sin^2 always regardless of the next symbol, but if the symbol is a zero the result will be zero throughout the duration of the zero. So in this case with this particular waveform, the way the lower frequency is constructed actually helps for a cleaner demodulation, since it provides the nulled response if the delay is the duration of the zero symbol (Maybe the waveform designer had this in mind which is why it has the peculiar form for the lower frequency).

The "Square Wave" output would be further filtered through an "integrate and dump" process prior to making a final symbol decision: sum over 12 samples in the symbol duration and decide if the result if greater or less than a decision threshold. If greater a "0" was demodulated, and if less a "1" was demodulated. Given the atypical longer duration between the symbols, four "1"'s in a row is the presence of an actual "1" symbol.

This structure can then be readily used with standard timing recovery for recovering the symbol clock boundaries. If the symbol clock is synchronized (coherent) to the waveform sample clock, the timing recovery is quite trivial. Two options show below:

The strong double frequency component can be used for symbol timing recovery for the case of the symbols no being coherent to the sample clock. If they are coherent then symbol timing recovery is simply aligning the sample edge to a start of a symbol which can easily done with a simple threshold detection at the Demodulated FSK output as shown in the block diagram (for lower SNR conditions it would be more robust to average the estimated timing position even such a threshold approach was used, such that the timing position does not change abruptly based on the result of any one threshold detection). Once detected, a modulo 12 counter starts, and every time the counter rolls over the symbol is selected at the output. The more robust timing approach would compare the roll-over timing clock and the FSK output threshold detector impulses to establish a timing error (another multiplier would do this), which would then be accumulated in a timing loop to adjust the actual start time of the cyclical timing recovery counter (again only need to do that for a very robust solution in the presence of low SNR, for the waveform as shown all that would not be necessary).

See this link where I detail non-coherent FSK demodulation approaches that further how a delay and multiply forms a frequency discriminator:

https://electronics.stackexchange.com/posts/293723/edit

With this graphic in mind; a delay multiply is a frequency discriminator:

Optional coherent techniques would have a 3 dB SNR advantage but require more processing. This could be done by matched filter correlation to the two symbols (have replicas of the two symbols and multiply and accumulate to each one to decide which symbol is observed). The delay and multiply followed by delay and add is notably simpler.


Covariance matrix of an adaptive filter input I run many times in equations containing the trace of covariance matrix of an adaptive filter input. But it is not really clear what it is.

For example in this paper the input covariance matrix is

$$\textbf{R} = E\left(\textbf{u}_i^*\textbf{u}_i \right)$$

where:

$\textbf{u}_i$ is row vector of input (seems to be regression vector in given time)
$^*$ denotes Hermitian conjugation (complex conjugation for scalars)
$E$ is not explained, but I would expect it is the expected value

Questions:

What does it mean?
How does the matrix look like, what is the shape of it?
How do I estimate it from input data $\textbf{U}$ - matrix where one row represents one input vector?

The Covariance Matrix is commonly defined as

$$\mathbf Q = E\left[ (\mathbf x -\mathbf\mu_{x})(\mathbf x -\mathbf\mu_{x})^*\right]$$
with $\mu$ denoting the mean value, i.e. $\mu_{x}=E\left[\mathbf x\right]$, and $\mathbf x$ being column vectors. The fact that you define the covariance matrix as

$$\mathbf{R}_i = E\left[\textbf{u}_i^*\textbf{u}_i \right]$$

indicates that your $\mathbf u_i$ have zero mean and are row vectors; that might be a result of the signal model you use, but it's pretty non-standard.

Furthermore, the product of a column vector with its hermitian is inherently hermitian (symmetric if real).

Now, one way to estimate an expectation value from multiple observations. Wikipedia has an article on estimating covariance matrices, so here's just the gist: the most intuitive estimator is the sample covariance matrix; adapted to your formula, that'd be

$$\hat{\mathbf R} = \frac1{N-1}\sum\limits_{i=1}^N \mathbf u_i^*\mathbf u_i$$


Thoretically, what makes multi-pitch detection so difficult? For multi-pitch detection, I mean for example pressing several piano keys at the same time. For signal pitch detection (pressing one key at a time), after FFT, mostly the frequency with the largest magnitude is desired frequency.  However, for multi-pitch detection, let's just say press two piano keys, the two frequencies with the first two largest magnitudes are not the frequencies of keys I pressed.

First of all:

  mostly the frequency with the largest magnitude is desired frequency

This is a very naive and flawed approach, as what you are measuring here is very different from the perceptual notion of pitch. Not to mention that you hit the limited resolution of FFT...

Some of the challenges one has to face in multi-pitch detection:

Octaves are difficult to tell apart. Play two Cs one octave apart, and the harmonics of the highest note reinforces the harmonics of the lowest note, rather than create new peaks in the spectrum, to the point that one could hear one single "C" note with a different timbre. When mixing several notes together, some of their partials can reinforce or cancel each other.
Pitch detection can no longer be seen as a problem of period estimation. The most successful monophonic pitch detection algorithms are time-domain and based on metrics like the average magnitude (or square) difference function, and the underlying assumption is that the signal has a salient periodicity - it is almost similar to itself when translated in time. These methods cannot be applied to polyphonic signals, as the period of a mixture of periodic signals is the LCD of the components' periods. Knowing this LCD does not give away the factors...
The fact that the number of notes is not known in advance increases the size of the search space, and makes it an ill-posed problem. You are no longer ranking hypothesis like "this is a C3" vs "this is a C#3", but instead, "this is a C3 and a G4" vs "this is just a C3" vs "this is a C2 and a C3". For example, take a single piano note. One could describe its spectra as a "single C3 piano note". Or one could describe its spectra as "a C3, a C4, a G4, etc... all played together by an instrument producing purely sinusoidal notes". Distinguishing among these two options requires prior knowledge of the spectral shape of a piano note.


Finding Repetetion rate of ns pulsed waveforms my application :
sensing a photo diode output and processing it further, the per-processed output would be a digitally converted pulse waveform with pulse width ranging from 1ns-150ns and pulse repetition frequency 1Hz to 50KHz,now i should find the repetition rate of the 1ns-150ns pulse

as sensing this kind of narrow pulse signals is difficult, my per-processing stage would expand input signal by a factor of 100 thus my signal would become a signal with pulse width 100ns-15000ns

to find pulse repetition frequency for which i thought FFT would do the job

problem statement:
my sampling rate would be 2MSPS - to acquire a 100ns pulse i require that much

so if i perform a 1024 FFT over 1024 samples i acquired i will not find the pulse repetition rate at all, in order to find the repetition rate i have to grab at least 1sec samples and do a very large FFT over it. which is practically foolish

so i am a bit confused in this situation how can i find the pulse repetition frequency of a 100ns pulse,through huge samples acquired, mean while i am of the idea to just using a pulse counter of micro controller.

problem statement: my sampling rate would be 2MSPS - to acquire a 100ns pulse i require that much

Nope. 100ns is 5 times as short as the interval between two sampling instances in your sampling device. If it does indeed measure the voltage accumulated over a capacitance, this might work – if that capacitance is as small as possible, which it often is, you might just miss pulses.

However, you're right, the signal you're interested in actually doesn't have a bandwidth of 10 MHz – the maximum frequency of interest is 50 kHz, so observing a system that indicates the presence of a pulse (not the pulse itself) with a 100 kS/s sampling rate should suffice.

If you know you just got presence information (i.e "there was an impulse during the last sample time, or not"), an FFT is a bit over the top. Simply count the samples between an observation.

The hard part actually would be the electronics to your photo diode – 1ns impulses imply that your system needs to have a bandwidth of 1 GHz – so you'd be deep in the world of RF/microwave circuitry. What you'll probably end up doing is amplifying the current coming from your diode, with something avalanche-based, or a sensitive microwave transistor. Don't forget to lay out your board following the needs of a microwave signal! You can't afford spurious interference here.

You'd essentially amplify the energy of the pulse. The voltage over a simple RC lowpass fed with that energy with a 90% discharge time of let's say 1/(100 kHz) would serve as a feasible indicator of activity; observe that, if you want with 2MS/s, and you'd have captured a signal from which you can directly see the impulses.

Notice that photodiodes with a > 1GHz bandwidths will be noisy, so you'll probably have some false impulses. Here, your FFT approach might actually tell the periodicities from the white noise.


Fourier transform 4 times = original function (2D and higher) The Signal Processing SE post linked below shows how the Fourier Transform applied 4 times to a 1D function returns the original function, i.e. F{ F{ F{ F{ g(x) } } } } = g(x)

Link to 1D case:
Fourier transform 4 times = original function (from Bracewell book)

Can someone prove that this does or doesn't apply to 2D or higher functions?
In particular, if I use this method on a grayscale image, will I get back the original image? I know edge effects could be problematic, so lets say the image is a simple repeating pattern so that it is periodic. And yes, numerical/computational FFT issues might mean the image is not exactly like the original, but even if they are similar that would be pretty cool to see....

The n-dimensional (discrete) Fourier transform is separable, that means the dimensions can be treated independently. So the property that the FT is a fourth root of the identity applies to n-dimensional transforms too. And the property is exact, so apart from a tiny numerical round off error you won't see any deviations.


Transfer Function definition To find the transfer function of a channel we say that it is
$$
H(s) = \frac{y(s)}{x(s)}|x(s)=0 for <0
$$
Why we do not define it like
$$
h(t) = \frac{y(t)}{x(t)}
$$

We do not define it in that manner. Note the difference between the time and Laplace domains.

Generally, we characterize a system by the response produced when the system is excited by an impulse response $\delta(t)$. Given the output $y(t)$ for such signal we ask which function $h(t)$ did we convolve with the $\delta(t)$ to produce it. Formulation:
$$\delta(t)*h(t)=y(t)$$

Here $*$ represent convolution and not a product. When we find $h(t)$ we can produce the output for any other input $x(t)$ by means of convolution:
$$x(t)*h(t)=y(t)$$

If we transform into the Laplace domain, a convolution is translated into multiplication:

$$\mathcal{F}\{x(t)*h(t)\}=\mathcal{F}\{y(t)\}$$
$$\mathcal{F}\{x(t)\}\cdot \mathcal{F}\{h(t)\}=\mathcal{F}\{y(t)\}$$
$$X(s)\cdot H(s)=Y(s)$$
$$H(s)=\frac{Y(s)}{X(s)}$$

It is a convention that the capital letters represent the transformed signals.  Also, $t$ has changed into $s$ because we are now in the Laplace domain. Since the convolution transformed into a product we may now divide both sides by $X(s)$. We could not do this on the time domain because there we had a convolution operator and not a product.


digital phase lead algorithm I have two sine signals A and B which are simultaneously sampled at 2500 Hz using ADC. Through oscilloscope, I find that signal B lags A by about 5 degrees. I would like to digitally lead signal B to sync it with A. How to do that in MCU?

Do they lag 5 degrees for all frequencies? Or only for a specific signal frequency?

If you want to shift a signal by 5 degrees for all frequencies you need some kind of all-pass filter similar to a Hilbert filter. However the Hilbert Filter whether FIR or IIR shifts a signal by 90 degrees not 5.

However, if you want to add or subtract 5 degrees at a specific frequency, you need to this

delay = phase/( 2*pi * f).

If your delay is an integer number of samples, you simply need to delay your leading signal by this number of samples. However, if the delay is fractionnal you need to delay the leading signal with a fractional-delay filter.

http://users.spa.aalto.fi/vpv/publications/vesan_vaitos/ch3_pt1_fir.pdf


As DCT Does Compaction of Energy, Can We Say IDCT Does Energy Spreading? Like energy compaction transforms, the inverse transforms like IDCT are energy spreading transforms?

For certain kinds of signals (voiced speech, most music, etc.), the DCT is likely to compact, and thus the IDCT to (re)spread that more compact DCT result; but that is not true in general.  An impulse in either domain will be spread in the other domain, and vice versa, for that spread result, when inverse transformed.


What does the one-tap equalizer used with OFDM accomplish? I've read that OFDM normally uses a one-tap equalizer.  It seems that a one-tap equalizer could only scale and delay a signal.  Why is this useful?

I thought multipath was handled through the cyclic-prefix and by choosing a symbol width that is large compared to the delay spread.  Why is an equalizer needed at all?

The OFDM signal as a whole is affected by frequency selective filtering. It is usually designed such that the subcarrier bandwidth is smaller than the channel coherence bandwidth. This yields you frequency flat fading for each subcarrier which can be described by a single complex multiplication and equalization can equally done with just one tap.

Equalization in OFDM usually happens in the frequency domain. Let $X$ be the vector of subcarriers in the frequency domain and $x=IDFT(X)$ the corresponding time domain signal. Then, if the impulse response of your channel is a vector $h$, the received signal $y$ will be $y=x*h$ where $*$ referes to convolution. Then $Y=DFT(Y)=X \cdot H$ with $H=DFT(h)$ and $\cdot$ refering to element-wise multiplication. However, this is only true if $y$ is the circular convolution of $x$ and $h$ while in reality the convolution is linear. Therefore, the cyclic prefix is prepended (it needs to be at least as long as $h$) so that the linear convolution contains the circular convolution.

Then, $Y=DFT(Y)=X \cdot H$ with $H=DFT(h)$ and $\cdot$ as stated before. Since each element of $Y$ and $X$ refers to one subcarrier it can be seen that each subcarrier is multiplied with just one complex element from $H$ and hence equalization just becomes $Z=Y/H$ where again $/$ referes to element-wise division. This describes zero-forcing equalization which might not be the best choice if additive noise and interferers are present but it transports the idea.

And this is nothing else but one-tap equalization for each individual subcarrier.


HVS models and rods and cones It is common to say that the YUV color model is inspired from the human visual system (HVS), knowing that we have many more rods than cones, and consequently, we are more sensitive to light than color. However, it does not seem to be that simple. Rods become active under low light conditions, and cones are not totally insensitive to light. These details seem to be ignored by the literature. My question is: is it reasonable to build HVS models for image processing only based on the number of rods and cones?

For instance, in the paper "A new color image quality measure based on YUV transformation and PSNR for human vision system" by Yalman and Ertürk (http://journals.tubitak.gov.tr/elektrik/issues/elk-13-21-2/elk-21-2-20-1111-11.pdf), a new quality metric called CQM is proposed. It consists in computing separated PSNRs on the Y, U and V channels, and combining them with a weighted average, whose weights are adjusted according to the number of rods and cones. But, once again, is it reasonable knowing that rods become active under low light conditions?

As far as i know the YUV color space was not inspired from HVS, but from technical compatibility reasons; to transmit new color video television broadcast signals over the existing black and white infrastructure.

Essentially, the Y component of a color image (which is captured naturally in RGB system) carries the black and white or the brightness information of the image and is transmitted in such a way that old TV equipments (those that do not have color receiving capabilities) can still receive and decode the Y component as if it were the old black and white signal and still be able to display it without any hardware modification.

While at the same time, new (back then in 50s, 60s) color TVs were capable of also receiving transmitted U and V components and were able to display color pictures by properly decoding the YUV signal using some sort of comb filtering in NTSC and PAL systems.

Furthermore, Y carries more bandwidth and information than U and V. There's a separate color carrier towards the end of the channel bandwidth (before the audio carrier).

Because of certain industrial standardization concerns, YUV and its variants (YPbPr, YCbCr) achieved wide range of acceptance and being used in almost every image and video electronic equipments and codecs (in its digital form).

It's true that Human sensititivity to color and brightness are different due to cone and rod structures. And researchers can take advantage of this information to optimize their metrics with respect to human visual system. The paper most probably assumes such a correlation.


AWGN and Rayleigh fading in slow flat fading channels in MATLAB I wanted to add Additive White Gaussian Noise (AWGN) and involve effect of slow flat fading on my transmitted signal in Matlab.

I created noise with zero mean and variance of one by

awgn_noise = random('norm',0,1,1,NoOfSamples);

Then I created Rayleigh fading by

rayleigh = sqrt( (random('norm',0,1))^2 + (random('norm',0,1))^2 );

As far as I understand, the rayleigh values will attenuate the whole signal while gaussian noise values will be added to each sample of transmitted signal as:

received_signal = (transmitted_signal*rayleigh) + awgn_noise;

Am I doing it in a right way?

As I described in this answer, it is generally assumed in the Rayleigh fading model that the energy is conserved on average (i.e. that the fading coefficient $R$ is such that $E\left\{R^2\right\} = 1$). As a result you would need to generate the Rayleigh distribution from two normal distribution with standard deviation of $1/\sqrt{2}$ (instead of with unit standard deviation). This can be achieved with:

rayleigh = sqrt( (random('norm',0,1/sqrt(2)))^2 + (random('norm',0,1/sqrt(2)))^2 );

or alternatively, since you are already using random from the statistics toolbox, with:

rayleigh = random('Rayleigh',1/sqrt(2));

Other than that, the expressions you have properly represent a slow flat Rayleigh fading channel for the transmission of a one-dimensional signal (such as BPSK, M-PAM, ...)


Is the book wrong? -- How to sketch frequency response obtained from H(z)? This question is directly related to digital signal processing, so I am asking here.
How to sketch frequency response obtained from $H(z)$?
I'm adding an example and its solution below. I did not understand some of the things for question 7.4, part D.
Any help is appreciated!

There are three points in graph. How and why the second from left is placed there? What's its purpose?
I have no idea how phase graph is created. What is the relation between x and y axis? What do they actually represent?
Why the book says "give me function of frequency" but the solution has $\omega/2\pi$ in x-axis?

Question 7.4:

An LTI system is described by the difference equation
$$y[n] = \frac{1}{3}(x[n]+x[n-1]+x[n-2]) $$
a. Determine the system function $H(z)$ for this system.
b. Plot
the poles and zeros of $H(z)$ in the z-plane.
c. From $H(z)$ obtain
an expression for $H(e^{j\hat{\omega}})$, the frequency response of
this system
d. Sketch the frequency response (magnitude and phase)
as function of frequency for $-\pi \le \hat{\omega} \le\pi$.
e.
What is the output if the input is
$$x[n]=4+\cos[0.25\pi(n-1)]-3\cos[(2\pi/3)n] $$

Solution:

As far as I can tell, the book is not wrong. (By the way, you're studying a 3-point moving average filter here, and it's useful to know about such filters.)

Your 1st question: That center dot on the Mag. and Phase curves serves no purpose that I can see. I suggest you not worry about that center dot.

Your 2nd question: The Phase curve is generated by evaluating the angle of the  $e^{-j\hat{\omega}}$ phase term, in the last line of the $H(\hat{\omega})$ expression in Part (c) of the Solution, as $\hat{\omega}$ goes from $-\pi$ to +$\pi$. As $\hat{\omega}$ goes from $-\pi$ to +$\pi$, the Solution's $\hat{\omega}/2\pi$ Phase curve freq axis label goes from -0.5 to +0.5.

(Note: as you study the Solution know that the $H(z)$ in Part (c) is not the $H(z)$ in Part (b). The $H(z)$ in Part (c) is the $H(z)$ from the Solution's Part (a).)

Your 3rd question: It's a good question. In the literature of DSP there are a number of different notations (different symbols with different dimensions) used to represent frequency. And those various notations certainly can be confusing at first. mmswe, have a look at the following web page, it might help you with your 3rd question:

http://www.dsprelated.com/showarticle/199.php


Digital Lowpass Filters, basic questions Forgive the newbie questions, I'm just starting out and wanted some basic information.

I have a signal coming from a digital gyro, most of the useful information is below 90hz but I am seeing significant noise spikes from mechanical motor vibration at 200-400hz (the spike moves depending on motor RPM).  This needs to be filtered in realtime to reject the noise with as little delay as possible.

When the sampling rate is increased, say from 1000hz up to 8000hz, does this increase the delay?
When sample rate increases, is the number of samples/cpu load per timestep for the same filter increased?
In the above situation, what is a good way to increase rejection of signals above 90hz while keeping the delay low?  ie: filter type, raise/lower sampling rate...etc.
An existing system I've seen uses an LPF followed by 2 bandstop/notch filters positioned over the noise spikes.  The designers say this 3 filter approach creates less latency than doing a single more aggressive LPF, does this sound right?

If you want as little delay as possible, I recommend notch filters at your resonance frequencies. If your resonance frequencies are fixed, it should be easy. Otherwise, a low-pass filter at about 100 Hz should do the trick.
No your delay depends mostly on the cut-off frequency, the order of the filter and the filter type. An order-2 IIR filter will have a smaller delay than a 30-tap FIR filter . Your sample rate will not affect (or barely) the delay.
Typically, the lower the cut-off frequency, the higher the delay and vice versa. Also, an order-2 low-pass IIR will typically have a smaller delay than an order-4 low-pass IIR for the same cut-off frequency.
If you can live with non-constant group delay, an IIR should do the trick.
Yes it does, but it assumes that you know what your resonance frequency is. You can maybe adapt your notch filter to the resonance frequency in real-time.. it is feasible but can be a lot of work.


How to Remove the Periodic Oscillations from a Signal The task that I have is to remove the annual and semiannual oscillation from a set of temperature measurements, taken over several years, by means of least squares method.

I found the method described here and since I have to wait for my colleague to first make some analysis on the data and pass them on to me, I decided to make a test run in Matlab, using a signal with two simple sine and cosine functions and random noise.

%% Creating the data

g=2.5;
t=linspace(0,g,g*365*24*60);              %n-years of 1min data
o=rand(size(t));                         %random noise
p=2*sin(2*pi*t)+1.5*cos(pi*t);           %we create two sinusoids
x=p+o;                                     %and add noise to them

%% The method

N=length(x);            %no of data
n=(1:N);
T=g*365*24*60;          %time that the data covers
f1=1/(365*24*60);       %frequency of 1 yr oscillations
f2=1/((365*24*60)/2);   %frequency of 1/2 year oscillations

a1=f1*T;
a2=f2*T;

A1=(2/N)*sum(x.*cos((2*pi*a1*n)/N));
A2=(2/N)*sum(x.*cos((2*pi*a2*n)/N));
B1=(2/N)*sum(x.*sin((2*pi*a1*n)/N));
B2=(2/N)*sum(x.*sin((2*pi*a2*n)/N));

C1=sqrt(A1^2+B1^2);                     %amplitude
C2=sqrt(A2^2+B2^2);

fi1=atan(B1/A1);                        %phase shift
fi2=atan(B2/A2);

v=(C1*cos(2*pi*t-fi1))+(C2*cos(2*pi*t-fi2)); %the reconstructed set of data
                                             for these two frequencies

r=x-v;    %I wasn't sure what was meant by removing, so i just subtracted

%% Visualization

figure                                  %now lets plot the data
plot(t,x)
xlabel('t[min]')
ylabel('Arb. units :)')
title('Data')

figure                           %and, original signal and noise separately
plot(t,p,'b')
hold on
plot(t,o,'-r')
xlabel('t[min]')
ylabel('Arb. units :)')
title('Data separated')

figure                    %and once again data (blue), reconstructed signal
plot(t,x,'b-')            %(green), and the final result afer subtraction
hold on                   %(red)
plot(t,v,'g--')
plot(t,r,'r:')
xlabel('t[min]')
ylabel('w/e')

The problem now is, that the reconstructed signal (v) does not resemble the original signal so much. Even when I remove the noise and work only with the pure sinusoid signal, what I get is not really a good representation of it.
So, what I'm most interested is:

Is this at all the correct use of this method? Am I using the correct equations?
Is there an error somewhere in the code? By error I mostly consider the wrong definition, or calculation of some of the variables.
In this "ideal case" (sinusoid without noise) should the method be able to reproduce the original curve (almost) exactly?

I hope I've laid out my problem clearly and understandably, if anyone has additional questions, please ask.
I know this is rather basic stuff, but my knowledge and experience in this area is very limited (up to non-exsistant), so please bear with me.
Thank you in advance for any help.

Their are some problems in the code and same algorithmic issues:
Code:

You are estimating the wrong frequency: f2=1/((365*24*60)/2); should be f2=1/((365*24*60))/2;

Algorithm:
Basically you are using a fourier series (http://en.wikipedia.org/wiki/Fourier_series) with those terms corresponding to the known frequencies. But you should be aware of the limitations of that approach:

You have to estimate the dc part (frequency 0), because your signals' dc value is not zero!
In theory this approach works only for full periods of each frequency. That's why f1 is estimated quite good, but f2 is not. To work arround that you need either more data, because the more periods are in the data the smaller will be the estimation error. And you should use a windowing function, to get rid of some side effects by cutting the signal at the beginnig and the end.

The improved code could look like this:

clear
doWindow = true;
scaleFreq = 1;
%% Creating the data

g=2.5;
t=linspace(0,g,g*365*24*60);              %n-years of 1min data
o=rand(size(t));                         %random noise
%we create two sinusoids
p1=2*sin(2*pi*t*scaleFreq);
p2=1.5*cos(1*pi*t*scaleFreq);
p=p1+p2;
x=p+o;                                     %and add noise to them

%% The method

N=length(x);            %no of data
n=(1:N);
T=g*365*24*60;          %time that the data covers
f1=scaleFreq*1/(365*24*60);       %frequency of 1 yr oscillations
f2=scaleFreq*1/((365*24*60))/2;   %frequency of 1/2 year oscillations

a1=f1*T;
a2=f2*T;

if(doWindow)
    window = hanning(numel(x)).';
else
    window = ones(size(x));
end
xWin = x.*window;%windowing
xMean = mean(xWin);%calculate dc
xWin = xWin-xMean;%subtract dc

buildFunkCos1 = cos((2*pi*a1*n)/N);
buildFunkCos2 = cos((2*pi*a2*n)/N);
buildFunkSin1 = sin((2*pi*a1*n)/N);
buildFunkSin2 = sin((2*pi*a2*n)/N);

A1=(1/N)*sum(xWin.*buildFunkCos1)
A2=(1/N)*sum(xWin.*buildFunkCos2)
B1=(1/N)*sum(xWin.*buildFunkSin1)
B2=(1/N)*sum(xWin.*buildFunkSin2)

if(doWindow)%do correction due to windowing
    A1 = 2*A1;
    A2 = 2*A2;
    B1 = 2*B1;
    B2 = 2*B2;
end

C1=2*sqrt(A1^2+B1^2)%amplitude
C2=2*sqrt(A2^2+B2^2)

fi1=atan(B1/A1)%phase shift
if(A1<0)
    fi1 = fi1 + pi;
end
fi2=atan(B2/A2)
if(A2<0)
    fi2 = fi2 + pi;
end

%the reconstructed set of data for these two frequencies
v1=C1*cos(2*pi*t*scaleFreq-fi1);
v2=(C2*cos(1*pi*t*scaleFreq-fi2));
v=v1+v2+xMean;%sum estimates (freq 1,freq 2, dc)
r=x-v;    %I wasn't sure what was meant by removing, so i just subtracted

%% Visualization
figure(1)                                  %now lets plot the data
plot(t,[p1;p2;p]);
legend('p1','p2','p');
title('single sine signals and sum');

figure(2)                           %and, original signal and noise separately
hold off
plot(t,[v1;v2;v]);
legend('v1','v2','v');
title('estimated sine signals and sum');

figure(3)                    %and once again data (blue), reconstructed signal
hold off
plot(t,x,'b-')            %(green), and the final result afer subtraction
hold on                   %(red)
plot(t,v1+v2,'g--')
plot(t,r,'r:')
xlabel('t[min]')
ylabel('w/e')
legend('measured','reconstructed','result');

figure(4)
hold off
plot(t,[buildFunkCos1;buildFunkCos2;buildFunkSin1;buildFunkSin2]),
legend('buildFunkCos1','buildFunkCos2','buildFunkSin1','buildFunkSin2');

Here you can switch the windowing and also test what happens if more periods are in the signal. Try, scaleFrequ=10 und doWindow=true and the result should be what you are looking for.


How to turn Morlet wavelet into second generation wavelet? I have a very specific application which requires to use the complex Morlet wavelet (or Gabor wavelet, if that's the name you use). I currently use the convolution theorem to compute wavelets transforms, but I'm reading everywhere that direct wavelet transforms are more computer efficient, because of the lifting scheme.

On the other hand, I searched for lifting scheme, and it's so much of a different paradigm that I have no idea how to use it with my wavelet, or whether it's even possible.

My question is: is it possible to use the lifting scheme with the Morlet wavelet in 1D and 2D? If so, can you recommend some litterature explaining the algorithm?

The Lifting scheme is special algorithm that allows a very efficient computation of the discrete wavelet transform (DWT), via special filters (QMF) and down/up-samplers. In order to compute a DWT the wavelet has to fulfill orthogonality constraints under translation and dilation.

The Morlet wavelet does not fulfill these contraints. Thus, you cannot compute a DWT with it und consequently you cannot use the lifting scheme. The Morlet wavelet is only applicable for computing a continuous wavelet transform. If you need a fast approximation you can have a look here: Fast Gabor-Like Windowed Fourier and Continuous Wavelet Transforms


What level of correlation coefficient suggests correlation? In Python (and Numpy), I recently did:-

a = np.random.randint(0,256,1000)

b = np.random.randint(0,256,1000)

and when I calculate the correlation coefficient $r$ using:-

np.corrcoef( a,b )

I get a value of 0.07944437.  Rounded, that's 0.08, 8% or about $1 \over 13$ th.

$r$ is said to typically range -1.0 thru 0.0 to +1.0.  I've never had a 0.0 and don't think that it's possible as random correlations can always be found between any two series. It's like asking "How small is small?" 0.08 doesn't seem that small. Yet clearly $a$ is not correlated to $b$ as both originate from the good quality Mersenne Twister deep inside Python.

So if $r = 0$ is empirically impossible, when does it actually become significant and suggest a real correlation?

Note 1.  There is of course the 68–95–99.7 rule. This covers a $3 \delta$ range of certainty. The atomic guys use $5 \delta$ to confirm a discovery. Do  signal guys do anything similar, and if so, how would they apply it?

Note 2. With the greatest of respect, How do I know quantitatively if the correlation of two time series is significant? seems ambiguous.

One way to look at it is that $r^2$ (yes, I do mean $r^2$, not $r$ or $|r|$) is the fraction of the variance of the target (the series A) that is explained by the estimator (a linear function $\alpha + \beta B$ of the series B).  Note that if Series A had variance $\sigma^2$, then the variance of $(A - (\alpha + \beta B))^2$ (what B cannot explain about A) is $(1-r^2)\sigma^2$.  So, if $|r| < \frac{1}{\sqrt{2}} \approx 0.707\ldots$, less than half the variance has been explained, and so one viewpoint is that a correlation of less than $70\%$ or so suggests that the linear estimator being considered here is not all that great.


Phase differences between DFT and FFT (in MatLab) I have implemented a (direct) DFT in MatLab (following script) and compared it to the built-in FFT routine. The magnitude response seems to be identical (excluding some possibly round-off errors), but in the phase I get different results.

I can't find the source of these differences and I would like some insight on that.

script (only the case of even vector length is presented here):

%% Calculate needed variables
DFT = zeros(length(data), 1); % Pre-allocate result vector
n = (1:length(data))/length(data); % Calculate the n/N factor
n = -1j * n * 2 * pi; % Calculate the (-j * n * 2 * pi)/N factor

dftLength = (length(data)/2) + 1; % Calculate the length of unique Fourier coefficients

% Calculate the coefficients
for i = 1:dftLength
    exponent = exp((i - 1) * n); % Calculate the exponent
    DFT(i) = exponent * data; % Multiply the signal with the exponent

    % Add the conjugate symmetric part of the spectrum
    if i ~= 1
        % Skip first bin (it's DC)
        DFT(2 * dftLength - i) = conj(DFT(i)); % Add conjugate symmetric bin
    end
end

You can see the resulting plot here

you should replace the line

n = (1:length(data))/length(data); % Calculate the n/N factor

with the line

n = (0:length(data)-1)/length(data); % Calculate the n/N factor

Note that in the standard DFT (and FFT), the indices for both $x[n]$, $X[k]$ range in $0 \leq n,k \leq N-1$.


Are IQ-samples independent from the carrier frequency? Since i am not sure if i have understood IQ samples correctly, i want to ask the following question:

If i look at inphase and quadrature in an IQ diagramme, the length of the arrow describes the amplitude of the signal, while the angle between the arrow and the x-axis describes the phase shift. That however would mean, that these samples could still be transmitted on different frequencys, maintaining the same amplitude and phase shift.

Is this right or am i misunderstanding something?

I'm not sure I understand your question, but here are some thoughts.

Say you have a passband (that is, with a carrier) quadrature signal: $$s(t) = x(t)\cos(2\pi f_c t) + y(t)\sin(2\pi f_c t).$$ Here $x(t)$ and $y(t)$ have bandwidth $B<<f_c$. If you sample this signal, the specific samples you get will obviously depend on $f_c$.

One almost never does that; for many reasons, it's better to remove the carrier. So, you can obtain the low-pass equivalent (complex envelope) signal $$s_{BB}(t) = x(t) + j y(t),$$ and sample that instead. Then, the samples will not depend on $f_c$.


why do we use $X(e^{j\omega})$ instead of $X(j\omega) $ in Discrete Time FT I am studying DT-FT. But I cannot figure out why we use $X(e^{j\omega})$ instead of $ X(j\omega) $ in DT FT

Thanks in advance..

The argument $e^{j\omega}$ emphasizes the $2\pi$-periodicity of the discrete-time Fourier transform (DTFT) of a sequence. Furthermore, if the $\mathcal{Z}$-transform $X(z)$ of the sequence $x[n]$ exists and if the unit circle $|z|=1$ is inside the region of convergence, then the DTFT of the sequence is simply obtained by evaluating $X(z)$ on the unit circle $z=e^{j\omega}$.

That's the same thing as with the continuous Fourier transform and the Laplace transform, referred to in your previous question.


Finding number of independent samples using autocorrelation I have a pressure signal (y) with 512000 samples and with a sampling frequency 5000 sample/sec. I am trying to find the number of statistically independent samples using autocorrelation. I used

autocorr(y)

I do not know if this is correct or how to get the number of independent samples

I found this question, but there is no valid answer

https://stats.stackexchange.com/questions/149517/autocorrelation-and-statistically-independent-samples

Thanks

The autocorrelation you obtained shows that the signal becomes decorrelated with itself for specific time delays, which are the zero-crossings of the autocorrelation. The autocorrelation is not zero except at these time instants.

You can also define a threshold $\tau$ and say that, if the autocorrelation is less than $\tau$, then for your purposes the signal is decorrelated for any time shift $\Delta$ for which the autocorrelation is less than $\tau$.

When calculating the autocorrelation in Matlab (or similar), the lag is calculated in number of samples. For example, in your case it seems like the autocorrelation is zero for a lag of 200. Since your sampling interval is $T_s = 1/5000 = 200\,\mu$s, then your signal becomes decorrelated for a shift of 200 samples, or a time shift $\Delta = 200 T_s = 40\,$ms.


How to recover frequency and damping coefficient from superposition of damped oscillators? I have a signal which is (more or less) the superposition of multiple damped oscillators.  There are relatively few oscillators involved, 2-3 would account for most of the signal.  They have periods which are not multiples of each other (ie not harmonics), and they all have different damping coefficients / time constants.  They also are not necessarily in phase initially.

I would like to recover the periods, damping coefficients and phases of the oscillators from the signal.

A synthetic example of a signal like this can be produced by the following python code:

from pylab import *
x = linspace(0,1000,num=1000)
y1 = 1 * sin(2*pi*x/11 + 0) * exp(-x/300)
y2 = 0.5 * sin(2*pi*x/17 + 1.2) * exp(-x/300)
y3 = 2 * sin(2*pi*x/50 - 0.25) * exp(-x/50)
y = y1 + y2 + y3
figure(figsize=(20,10))
plot(y)

This looks pretty similar to the real signal qualitatively, including having some beat frequencies, and a pretty large initial low-frequency component.

This seems like a pretty common problem when analyzing a physical system with multiple coupled oscillators, eg a vibrating object with multiple resonances.

UPDATE Using direct curve-fitting works pretty well, even for strongly damped signals - depending on the choice of initial conditions, that is:

from scipy.optimize import curve_fit

def func(x, a, b, c, d):
     return a * sin(2*pi*x/b - c) * exp(-x/d)

p, pcov = curve_fit(func, x, y, p0=(1,100,0,1000))
print(p)

# [  1.78453351  51.53572084   0.05566558  55.10643076]

This recovers the largest-amplitude (but very strongly damped: time constant = 1 period) part of the signal.

If you are absolutely sure that you have a signal that is the superposition of damped oscillators then you can simply track their evolution in time.

Fs = 2000; % Sampling frequency (Hz)
T = 4;     % Time duration (seconds)
t = 0:(1./Fs):(T-(1./Fs)); % Time vector (in seconds)
p = 2.0.*pi.*t; % Phase vector
S = 1.0*sin(440.*p).*exp(-t) + 0.25 .* sin(800.*p).*exp(-t*10); % Signal

This signal looks like:

And via a Short Time Fourier Transform (STFT) (surf(stft(S))) it looks like this:

Granted, if you have beating frequencies, the task is getting a bit more challenging. For instance:

S = 1.0*sin(440.*p).*exp(-t) + 0.25 .* sin(442.*p).*exp(-t*2) + sin(800.*p).*exp(-t*10);

Would show up as:

But once you see oscillations on the damped decay, you know you have beating and you can recover how close the beating pair is by measuring the beating frequency and how quickly it fades (to derive its decay).

In a similar way (i.e. using the STFT) you can also recover the phases of the sinusoids.

Hope this helps.

UPDATE:

The STFT of the signal that was shared via pastebin is capable of distinguishing three components in this signal. This is as follows:

%Given a signal in s and that the signal package is installed (in Octave)
[S,F,T] = specgram(s);
surf(abs(S));

This returns:


Concepts: Ideal and Non Ideal Channel Resulting in Inter Symbol Interference I am following the book: Blind equalization and Identification by Zhi Ding and Ye Li. In Chapter 2, the concept of T spaced equalizers is presented. It is mentioned that the output of the channel is $x[k] = \sum_{i=-\infty}^\infty h[i]s[k-i] + w[k]$ where $h$ is the impulse response, $s$ are the symbols / source input and $w[k]$ is iid noise. Then it is defined that *

  "When the channel is non-deal, its impulse response $h[k]$ has more
  than one nonzero component. Consequently, undesirable signal
  distrotion as the channel output $x[k]$ depends on multiple symbols in
  $\{s[k]\}$. This phenomenon is known as inter-symbol interference
  (ISI)."

*

I do not understand why the impulse response has to be zero when it is ideal channel - no ISI. If $h$ is zero then $h[i].s[k-i] = 0$ where $dot$ is the multiplication operator. Or am I mistaken and in the equation $h$ is convolved with $s$ and not multiplied and essentially impulse response will be zero for ideal channel. What is the concept and reason? I do know that convolution in time domain is equivalent to multiplication in frequency domain. But the above representation is in time domain. A brief explanation will work and then I can follow on from the key points mentioned in the answer. Thank you.

The book doesn't say that the impulse response must be zero for an ideal channel. It says that an ideal channel has exactly one, and not more than one, non-zero component, i.e. the ideal channel's impulse response is an impulse, which means that the signal is only delayed but not distorted.


Calculating an exponential curve I'm developing a synthesiser and have finished the LFO's.   I'm wanting to provide a knob on the UI for setting the rate.   The minimum value of the LFO will be .01 Hz and the maximum value will be 3000 hz.   I'll have 100 increments in the knob and I want to apply an exponential curve to the values (Lfo rate).

What is the formula/function for calculating this?

Set $\Delta=\frac{1}{99}\log_{10}\left(\frac{3000}{0.01}\right)$

Then the $n_{th}$ knob position should be

$p_n=10^{n\Delta-2}$

Where $n$ runs between 0 and 99 in integer step.


How can an audio segment report an integrated loudness (LUFS) greater than 0 According to my very rudimentary understanding of the ITU-R BS.1770 standard, I thought that the LUFS (LKFS) value for any audio signals should not ever go positive.
However, there is a 5-second segment from an audio file I have whose integrated loudness level is +1.28 when measured using ffmpeg, and +1.44 when measured using pyloudnorm:
(stereo, source audio resampled at 48000 Hz, 16-bit signed little-endian integer waveform PCM)
[Parsed_loudnorm_0 @ 0x7fc245008680]
{
    "input_i" : "1.28",
    "input_tp" : "0.86",
    "input_lra" : "2.90",
    "input_thresh" : "-10.16",
    ...

How is this possible? Could this be an issue related to the implementations (i.e. ffmpeg and pyloudnorm) and not the algorithm itself?
I concede that I don't fully understand the math in the algorithm, and so if the algorithm does allow a positive value in certain situations, please explain.

I did download this annoying video and analyzed the part with a vuvuzela. You should take a note that:

This part is clipped
The audio is stereo

And that pretty much solves the mystery.

If you calculate the loudness for each of the channels separately you will get approximately -1.4 LUFS. Since this is a stereo track you add those levels together. For incoherent sources, you would roughly expect a 3dB of boost. Lo and behold, the total loudness calculated from the stereo track is approximately 1.6 LUFS (+3.2).
If you think about it, this makes perfect sense. What happens if you take one of those clipped channels and play it back through a speaker? You gonna get some SPL. Now, add a second speaker and play the other channel, you would expect to be even louder (3dB louder to be precise). Therefore LUFS accounts for that.


nature of output of the addition of white noise with colord noise I want to know that the addition of a uncorrelated time series x with a correlated time series, say y generates a time series but is that correlated or uncorrelated output? For example: x is White Gaussian noise and y is pink noise, then z = x+y Looking at the plot which is the output of the noisy time series, it is hard to say whether z is correlated or uncorrelated. Is there a rule about the nature of the output?

Are $x$ and $y$ uncorrelated? In this case the autocorrelation function of $x+y$ is just the sum of the autocorrelation functions of $x$ and $y$, which if one of them is not white will lead to its sum also not being white.
In general, for stationary processes you have $\phi_{zz}(\tau) = \phi_{xx}(\tau) + \phi_{yy}(\tau) + \phi_{xy}(\tau) + \phi_{yx}(\tau)$ where $\phi_{xy}(\tau) = \mathbb{E}\{x(t)y(t+\tau)\}$. Moreover, since the process is stationary you also have $\phi_{xy}(\tau) = \phi_{xy}(-\tau) = \phi_{yx}(\tau)$. So in general, the autocorrelation of $z$ depends not only on the autocorrelations of $x$ and $y$ but also their cross correlations.


What is the optimal way to select by color? I need to select all pixels in a image that are similar in color to a given one.

My current approach is to get the hsv values for this pixel and calculate the distance to each pixel in the image being

distance = abs(ref.h-pix.h) + abs(ref.s-pix.s) + abs(ref.v-pix.v);

And select the pixels bellow certain distance.

Is there any other preferred method?

The Wikipedia page Color difference gives several formulas, the most advanced being CIEDE2000. They are all based on Euclidean distance $\Delta E$ (straight line distance) in the approximately perceptually uniform CIE Lab color space, with various corrections.


High-pass filtering pre-processing before computing audio features I want to compute acoustic features on a set of audio files. These features can be classical audio features as spectral flux, but also acoustic indices as proposed in the seewave R package ).

I need, in a pre-processing step, to remove low frequency noise (approximatively inferior at 200 Hz), that could be caused for example by wind or a distant road, and are not appropriated for the analysis.

As my features can be computed both from temporal and spectral values extracted from the audio signal, I am considering the use of a temporal high pass filter. This approach seem to me more adequate than a filtering on the spectral domain (with a reconstruction of the signal by an inverse Fourier transform). Am I right?

Therefore I am looking for an high-pass digital filtering method that will remove low frequencies, but preserve as much as possible the high frequencies. If it is possible, I would like that the amplitudes of the fft bins above the cutoff frequency could stay the same after the filtering.

What is the best method for this application? Is it relevant to use a temporal filter rather than a filtering on the fft bins?

Concretely, I will use a python implementation (with scipy). I am now considering the use of a butterworth filter with the filtfilt function. Is this method appropriate? Is there a better solution for this issue?

Ant help or suggestion will be appreciate.

Edit

Sampling rate of the audio files: 48 kHz
Low frequency band to remove: around 0-300 Hz

I agree that a large-order (12 taps+) high-pass Butterworth filter will be sufficient.

Why?  Well, Butterworth filters are commonly used in professional audio-editing applications to do filtering, because their frequency response is designed to be as flat as possible in the passband. This means that the filtered audio is of higher quality/fidelity, as there are less frequency "ripples" or "artifacts" in the passband than other comparable filter types (Chebyshev, elliptic). In other words, Butterworth has a more linear phase response in the pass-band.

The downside is that the Butterworth filter has a slower roll-off, and will require a higher-order to implement a particular stopband specification. But this shouldn't be an issue for you.

To answer your second question, this is why the above solution is preferable to using an FFT and zeroing out the bins: Why is it a bad idea to filter by zeroing out FFT bins?