Estimating degree of motion in video frames I have a video (in webm format), and I want to produce a graph that shows the degree of motion in this video as a function of video time. I'm sure there are many different ways to define "degree of motion", but I envisaged something that could be based on the average difference between pixels in consecutive frames, or something like that (suggestions are welcome). Essentially what I need is a graph whose magnitude is large in cases where there is a very big change in the frames and small when most of the frame stays the same except for a small section.

What tools might I use to achieve this? Thanks in advance.

Probably you could try optical flow with some key-points not pixels if your purpose is to estimate motion. I thing meaningful thing would count average values of keypoints, overcomes some threshhold or filter (to eliminate barelly moving keypoints). However there would be some problems with having constant amount of meaningfull keypoints (keypoints tend to dissapear over time and very little chance that same amount newly generated keypoints are on moving objects). If that is what you want. You can estimate amplitudes of moving keypoints.
Tools could be any CV library like OpenCV. Pixels are not easily trackable without being extracted as keypoints.

On the other hand if you want to estimate how much pixels moving on image and estimate movement results from codec movement compression you could extract frame by frame and ignore I type frames and count how much pixels are on B and P type frames (I am recalling from h264 not sure how these frames are expressed in vp8), which are just moving picture parts from previous/next frame. You then can have estimate how much pixels are moving. Your call how you will implement this.


Why would the audio from .mp4 files be audible on some devices but not others? I don't know if this question is best suited for Signal Processing, Sound Design, or Video Production. If someone knows for sure that can do so, please migrate this question to the best forum for it (I'm guessing "Signal Processing" for now).

I have posted some mp4 files on the Internet, both on YouTube and on a dedicated facebook page. Most people have no problem with seeing the video and hearing the audio emanating from these files.

Two people, though, have complained (or mentioned, anyway) that they cannot hear the audio. One of them was using a Windows desktop, and the other an iPhone.

What are the reasons why some may not hear audio, while others can - when they can hear audio otherwise on their device?

The files in question (half a dozen of them) are here.

For a direct link to one particular one, you can try this.

Or, from YouTube, here's one.

Is it anything on my end that I can "fix," either in-place or locally and then re-upload? It'd be hard to test when I had accomplished my goal of making them audible to all, as they work for me already...

When such is the case, it's most probably a broken flash or incompatible browser for those who cannot hear the sound of your videos. However if those two can hear all other internet videos but only yours, then this is really a complex issue. Considering the fact that there are so many different platforms, browsers, devices, drivers and media formats it will be really hard to find out why...

for example I can not even run the videos you have uploaded on my windows 8.1 laptop with firefox browser. It's probably because I'm using a flash blocker to prevent unwanted streaming of flash videos whenever a page is loaded. Also firefox is known to casue erratic behaviour for flash videos through its latest versions...


Processing signals with phantom fundamentals Is there a formula to calculate the frequency of a phantom fundamental based on the harmonics, and would this formula be applied directly to the energy bins of an FFT?

Methods for estimating a fundamental pitch frequency (missing or not) from an FFT result containing a series of harmonics include finding maxima in a cepstrum, the harmonic product  spectrum method, and looking for clusters in histograms of the distances between frequency pairs of higher order frequency peaks that fit within psycho-acoustic filter bandwidths..  These are all more estimation algorithms that post-process an FFT result, rather than a simple formula.


Making a 3D point cloud from multiple RGB-D images I have multiple (4 kinect) cameras that give RGB-D (color and depth) information of the same scene from different points of view. I want to make a 3D point cloud out of these depth maps. I can get a 3D point cloud for each camera using a function available via the SDK. but my question is how can I merge these points in the same coordinate system to get a single point cloud? what are the available methods? are there any libraries available that perform this in python ?
My best guess is to use one of the cameras' coordinate system as reference and transform every point into that coordinate system. but how can I figure out the representation of one camera's coordinate system in another camera's coordinates so I can perform the rotation and translation ?
I am new to computer vision so please give some references I can use for the methods you mention.thank you.

Simple answer: Use LiveScan3D. Should be easy enough.
Longer answer: You are looking for a way to 'register' multiple point clouds. This is an extremely well studied problem in computer vision. In your case I would consider two scenarios:

Target-based (photogrammetry-like): Place some targets in the world such as markers (see OpenCV Aruco markers). Estimate the marker pose in each camera image (the camera center is set to be the origin). Once the marker pose is computed in all views, you could transform each 3D cloud into the coordinate frame of one of the cameras (choice is arbitrary).

Markerless: The 3D geometry itself is usually rich enough to relate the views together. This is called registration. You are indeed looking into a multiview registration problem. The typical pipelines follow two main steps: (i) pairwise registration; you can use for example FGR for this. Then these individual pairwise transformation estimates have to be linked together. If your view-graph in which there exist an edge for each pair of cameras seeing common structure, is simple, you could compute the minimum spanning tree and chain the transformations along the tree to compute all the coordinate frames with respect to the first one (free to choose).

Note that, with the advances in deep learning we can now solve this problem in an end-to-end manner, efficiently:
Gojcic, Z., Zhou, C., Wegner, J. D., Guibas, L. J., & Birdal, T. (2020). Learning multiview 3D point cloud registration. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (pp. 1759-1769). https://arxiv.org/pdf/2001.05119.pdf
The first option requires you to physically place fiducial tags in the scene but is easier to implement. The second option is more automatic, but might be harder to get to work.


Wav To Spectrogram, Back To Wav Currently, I'm writing a Python script, which should do the following:

read an audio file respectively a wav file via scipy.io.wavfile.read().
calculate the spectrogram of given wav file.
write the data from spectrogram back into a wav file.

Here's a bit of Code:

# Define FFT params:-------------------------------------------------------
windowSize = 512
shiftSize = 160
nFFT = 1024
window_py = signal.hamming(windowSize)
nOverlap_py = windowSize-shiftSize

# Load wav file into memory:------------------------------------------------
fs_rate,s_orig = wav.read('demo.wav')

# Type Casting:-------------------------------------------------------------
s_orig_py = np.asarray(s_orig,dtype=np.float64)

# Spectrogram:--------------------------------------------------------------
Fpy,Tpy,Spy=signal.spectrogram(s_orig_py,fs=fs_rate,window=window_py,
noverlap=nOverlap_py,nfft=nFFT,detrend='constant',return_onesided=True,
scaling='spectrum',mode='complex')

#---------------------------------------------------------------------------
P_py = np.angle(Spy)                    # Phase extraction:
X = np.absolute(Spy)                    # Needed for neural network!
X1 = X*np.cos(P_py)+1j*X*np.sin(P_py)   # "orig." spectrum. Needed for resyn

# Resynthesize to wav:------------------------------------------------------
X1 = np.append(X1,np.conjugate(X1[-1:1:-1,:]),axis=0)
x_opt_py = synthSpectrogram(X1,shiftSize,nFFT,window_py,nOverlap_py)
wav.write('demo.wav',fs_rate,x_opt_py)

But there's is a huge problem: The data in "Spy" is not useable. When I try to write the data back in a wav file, the result is noisy respectively there's nothing at all. Furthermore, I have a Matlab file, which does the same as the code above and it works just fine and in both cases the parameters are the same.

The values in of S in Python and Matlab are not even close and I don't understand why. I get that they can't be identical due the fact that both functions use a different algorithm to compute the FFT but as I mentioned before they are not even close.

Matlab-Values

Python-Values

A spectrogram only contains magnitude info, not phase info.

That means half of the original's signal information is lost; thus, the original signal cannot be recovered.


What is the inverse Z transform of this: $X(z) = \displaystyle \frac{1}{z}{\left(1-\frac{z^2}{4}\right)\left(1+\frac{1}{z}\right)\left(1-z\right)}$

Using partial fractions expansion i came up to this:

$\displaystyle \frac{1}{X(z)} = \frac{-\frac{z}{3}}{1-\frac{z}{2}} + \frac{\frac{z}{3}}{1+\frac{z}{2}} + \frac{\frac{2z}{3}}{1-z}$

I tried but i can't go on any further.

You just need to multiply the terms and sort by powers of $z$. If you do so you should get

$$X(z)=\frac{z^2}{4}-\frac54+z^{-2}\tag{1}$$

So the work is in multiplying out the terms, not in actually figuring out the inverse $\mathcal{Z}$-transform, because from $(1)$ you can immediately write down $x[n]$:

$$x[n]=\frac14\delta[n+2]-\frac54\delta[n]+\delta[n-2]\tag{2}$$


A little bit better than zeroing bins in order to filter The topic of zeroing frequency bins before IFFT has been discussed here : Why is it a bad idea to filter by zeroing out FFT bins?
It was very helpful to understand why it is a bad idea.

What do you think about, instead of zeroing, multiplying by a window (in the frequency domain!) which is 0 at the centre and smoothly goes to 1 outside of a small interval ?

t = linspace(0, 1, 256, endpoint=False)
x = sin(2 * pi * 3 * t) + cos(2 * pi * 100 * t)
X = fft(x)
window = 1 - scipy.hamming(128)
X[64:192] *= window
y = ifft(X)

Then it would not be a "rectangular" zeroing, but a "smoothed zeroing" in the frequency domain. What do you think about such a filter ?

The Fourier transform of your window-shaped frequency response is its impulse response.  Evaluate the transform of your window function, look at its length above your desired noise floor, adjust the filter response as required to make the impulse length appropriately short, add zero-padding of this length to your original signal before the FFT, and you approach classic convolution filtering.

The added zero-padding is required to remove or reduce artifacts due to circular convolution using FFT/IFFT fast convolution.  Knowing the length of the impulse response is required to know how much zero padding is needed.  Smoother (with wider transition bands) frequency responses tend to produce shorter impulse responses.  And using smooth window shaped functions tends to imply smoother wider transition bands (although there are many other methods that may do better).


Can time-invariance be determined from a given a transfer function? I've the following function.
$$ G(z) = 2 + \frac{-1+5z^{-1}}{(1-0.5z^{-1})(1-z^{-1})}$$
Calculating it's inverse using $\mathcal Z$-Transform, I get the following function:
$$g[n] = 2\delta[n] + 8u[n] - 9(0.5)^{n}u[n]$$
where $u[n]$ is the unit step and $\delta[n]$ is the impulse.
I tried to do the usual approach to find if a system is time-invariant or not.
That is shifting
the input and compute the output $y_1$ and shifting the output $y_2$, compare both; if equal then
it's time-invariant, else it isn't time-invariant. However, I don't know what could be the input in this case.
So, I'm supposing that $g[n]$ is the impulse response of the system. But I couldn't find a way to tell if the system is time-invariant or not based on it's impulse response. Maybe I'm missing some
key idea here.

The question can only be answered if it is clarified what it is that the function $G(z)$ describes. If $G(z)$ is the system's transfer function then we're done immediately, because only linear time-invariant (LTI) systems can be characterized by a transfer function of that form.
If $G(z)$ is the transfer function, the output sequence $y[n]$ is given by the convolution of the input sequence $x[n]$ and the impulse response $g[n]$, which is the inverse $\mathcal{Z}$-transform of $G(z)$. Note that in order to compute the inverse $\mathcal{Z}$-transform you must know the region of convergence (ROC), or you must be given the additional information that the system is causal.
Having computed the impulse response $g[n]$, the output is given by the following convolution sum:
$$y[n]=\sum_{k=-\infty}^{\infty}g[k]x[n-k]\tag{1}$$
It's trivial to show that the response to a shifted input $x[n-n_0]$ is given by $y[n-n_0]$. Hence, any system whose input-output relation is described by a convolution $(1)$ is time-invariant (and, obviously, also linear).


Generating good synchronization sequences As part of my work I am putting together a simple bursty BPSK communication system.  I want to make the demodulation easy, so I am going to put a synchronization sequence at the beginning of each burst to help the receiver lock onto the frequency, phase, and timing.

In order to get good synchronization information and avoid false positives the autocorrelation of the synchronization sequence should be as close to a delta dirac as possible.  There may also be additional constraints, like starting the sequence with all 1's to help the receiver get frequency lock quickly.

How does one go about generating such a sequence?  Is there a systematic way of doing it, or do you use common sense and trial and error?

A synchronization sequence generally needs the property that its autocorrelation function resembles an impulse.  There are two possible autocorrelation functions that can be considered. For a (real-valued) sequence $x$ of length $N$, the periodic autocorrelation function is
$$R_x[n] = \sum_{k=0}^{N-1}x[k]x[k+n]$$
where the sequence is assumed to extend periodically so that
$x[k+n] = x[k+n\bmod N]$.  In this case, one commonly used solution to
choose $x$ to be a maximal-length linear binary feedback shift-register
(LFSR) sequence of
period $N = 2^m-1$ whose periodic autocorrelation function is
$$R_x[n] = \begin{cases} N, & n \equiv 0 \bmod N,\\ -1, & n \not\equiv 0 \bmod N,
\end{cases}$$
which resembles an impulse train of period $N$ except that the "out-of-phase"
value of the periodic autocorrelation function is not identically
zero, but is small (and constant)
compared to the "in-phase" value of $N$.

The other form of autocorrelation function is the aperiodic autocorrelation
function which is given by
$$C_x[n] = \begin{cases}\sum_{k=0}^{N-1-n}x[k]x[k+n], & 0 \leq n < N,\\
\sum_{k=0}^{N-1+n}x[k-n]x[k], & -N < n < 0,\\
0, & |n| \geq N.
\end{cases}$$
Now, for binary sequences where $x$ takes on values in $\{+1,-1\}$,
the value of $C_x[n]$ alternates between even and odd integers as
$n$ varies from $-(N-1)$ to $(N-1)$, and so a "flat" out-of-phase
aperiodic autocorrelation function is impossible to achieve.
The closest we can hope to come is to find sequences for which
$|C_x[n]| \leq 1$ for $1\leq n < N$.  Sequences with this property
are called Barker sequences, and the longest known Barker
sequence is of length $13$.  it is widely believed, though not
proven, that Barker sequences of longer lengths do not exist,
and numerical studies show that if such a sequence exists, it
it must be so long that it will not be be usable in any practical
scheme as a synchronization preamble.


How to quantize filter coefficients to 10 bits I am trying to learn digital filter design. One of the concept is quantization of filter coefficients. I used Designfilt to design an lowpass IIR filter and obtained coefficients. Now I want to quantise the filter coefficients to 10 bits. How can I do that.
filterObj = designfilt('lowpassiir', 'PassbandFrequency', .45, 'StopbandFrequency', .55, ...
'PassbandRipple', 1, 'StopbandAttenuation', 60);
fvtool(filterObj);

[b,a] = sos2tf(filterObj.Coefficients);;

This is the code that I used to get lowpass filter.
Here b, a are numerator and denominator coefficients.

Don't quantize the numerator and denominator coefficients of the total transfer function. Use the coefficients of the second-order sections - sos in Matlab - and quantize those. You will usually implement an IIR filter using second-order sections, especially in a fixed-point implementation.
The first thing to do is to scale the coefficients of the second-order sections (SOS). Note that you can scale each section individually. Scaling means multiplying the numerator coefficients of each SOS such that you minimize the chance of overflow in that section. The scaling depends on the chosen structure of the SOS, and on the filter coefficients. A detailed discussion on scaling can be found in the textbook Digital Filter Design by Parks and Burrus. If you're interested in the design and implementation of digital filter you should definitely read that book.
After having scaled the coefficients of each SOS, you quantize them by simply rounding them to the closest number on an equidistant grid consisting of $2^N$ numbers, where $N$ is the number of bits.
Example: The coefficients of one of the SOS are given by
b = [1, 1.621783996, 1]
a = [1, -4.030702997e-01, 2.332661953e-01]
You decide to scale the coefficients of that section by a factor of $\frac12$, which of course only affects the numerator coefficients:
bs = [0.5, 8.10891998e-01, 0.5]
You use $N=10$ bits to quantize the coefficients. Since the coefficients are in the range $[-1,1)$, we can multiply them by $2^{N-1}=2^9$, round the result to the nearest integer, and then multiply by $2^{-9}$ to scale back to the original range. In Matlab/Octave this is easily done as follows:

b = [1, 1.621783996, 1];
a = [1, -4.030702997e-01, 2.332661953e-01];
scale = 0.5;
b = b * scale;
N = 10;
range = 2^(N-1);
ac = round( a(2:3) * range ) / range;
bc = round( b * range ) / range;

The quantized coefficients are given by
ac = [-4.0234375e-01, 2.32421875e-01]
and
bc = [0.5, 8.10546875e-01, 0.5]
Or, in $10$ bit binary two's complement form
ac = [1100110010, 0001110111]
bc = [0100000000, 0110011111, 0100000000]


Digtial FIR Impulse response & transfer function I am currently working the figure through below.
as it is an FIR Filter i have worked out using convolution that the output is 4,2,4,6,0,0.
i am trying to obtain the 'z' domain transfer function of the filter below using h(n)
using the $h(n)$ 2,-1,3
I have started to obtain the $z$ domain transfer function with
$y(n)=2x(n)-1x(n-1)+3x(n-2)$
any help?

You can solve this with high-school algebra and the time delay property
$$ \mathcal{Z}(x(n-k)) = z^{-k} X(z) $$
Applying this to your expression
$$y(n)=2x(n)-1x(n-1)+3x(n-2)$$
we get
$$Y(z) = 2 X(z) - X(z)z^{-1} + 3X(z)z^{-2}$$
$$\frac{Y(z)}{X(x)} = 2z^{-2} - z^{-1} + 3$$
The denominator is null because the filter is not recursive


How to make the harmonic complex tones sound comfortable? I'm trying to create a complex tone, which is consisted by several sine waves with a 220Hz fundamental frequency and 6 multiplied harmonics (at 440Hz, 660Hz, 880Hz, 1100Hz, 1320Hz, 1540Hz).
My code below generated a tone. However, it sounds very uncomfortable, and its waveform and spectrogram picture also look much intensive.
fs = 44100;      % [Hz]
T = 0.5;          % [s]
x = 0:1/fs:T;   % [s]

f = 220;        % [Hz]
y = 1*sin(2*pi*f*x);
f1 = 0.2*sin(2*pi*440*x);
f2 = 0.2*sin(2*pi*660*x);
f3 = 0.2*sin(2*pi*880*x);
f4 = 0.2*sin(2*pi*1100*x);
f5 = 0.2*sin(2*pi*1320*x);
y1 = y + f1 + f2 + f3 + f4 + f5
fname = 'sin_220.wav';
audiowrite(fname,y1,fs);

I found that one research group have shared similar tone materials, in which, the sine waves were synthesized regularly and sound pleasant.

My question is how to make adjustments to the code so that it looks like this.

In Octave: plot(y1);

audiowrite will effectively clip floating point sample values to range -1 ... 1. This will sound unpleasant. You can divide your data by its maximum absolute value, before writing the audio file, to normalize it to that range.


How to compute the free distance in TCM (Trellis-Coded Modulation)? I would like to be able to compute the so-called "free distance" for a TCM with a rate 2/3 convolutional code (two inputs, three outputs) inside that I plan to use for 8PSK. My previous assumption was that I could simply compute distance between the "null sequence" and the next closest sequence because the code is linear. But that does not seem to work because this computation gives me a squared free distance of supposedly 5.17157 for an 8-state code that actually performs worse in AWGN channel simulations than the 8-state recursive Ungerboeck code for which I compute the squared free distance to be 4.58579 (which is correct in this case).

Any pointers on how to compute the free distance for TCM would be greatly appreciated.

The TCM for which my free distance calculation seems to fail is this one:

in1 ----------*---------------*-- out2
              |               |
              |               V
in0 -*--------|--------*---->(+)- out1      (+) is a binary plus (xor)
     |        |        |
     V        V        V                    [D] is a delay element/flip-flop
    (+)->[D]-(+)->[D]-(+)->[D]-*- out0
     A                 A       |            V and A are supposed to be
     |                 |       |            arrow heads
     `-----------------*-------´

symbol mapping:
binary "out2 out1 out0" -> exp(2i*pi*out/8)
where out_k is the k-th bit of out counting from zero as the LSB

The free distance of a Trellis Coded Modulation system implemented either on BPSK or QPSK is a straightforward calculation. You'd have to compute the free distance $d_{free}$ of the convolutional code which is known as the Minimum Squared Euclidean Distance (MSED). Now, if the code has set-partitioned labeling to achieve coding gains over conventional labeling, then the minimum distance is given as $$\text{min}(d_{free},d_{ic})$$
where $d_{ic}$ is the minimum squared Intra-coset distance which can be calculated from the mapping. Having said this, the MSED of the convolutional code can be computed the same way its $d_f$ (the hamming free distance) is computed. You need to do a trellis search and replace hamming weights with squared euclidean minimum distances with the all zero-sequence mapping (this is valid because any convolutional code is linear and if the mapping is BPSK or QPSK, every symbol has a decoding Voronoi region that has the same area and the same shape as every other sequence). I believe you've implemented your trellis search perfectly. Having said this, the Ungerboeck code cannot be performing better than a code which has a higher Euclidean $d_{min}$. Suggestions:  1. Either you're computing one of the distances wrongly2. Are you sure about the minimum distance of the Ungerboeck code?


Given a spectrogram, how do I calculate loudness over time? I want to calculate the loudness of a given sound over time (like, in a graph) from its spectrogram. Can I do this from the spectrogram or is it better to do from the raw sound?

Yes, in fact it's not a bad idea at all to use some kind of spectrogram representation to calculate a loudness envelope. Depending on what exactly you need you can either just sum the spectral magnitude over all frequencies at a given time, or first apply some kind of perception based weighting functions
to make your loudness measure more related to actual auditory perception.

So if you give some more detail about the application you have in mind you will likely receive a more accurate answer.


Band-pass filter doubles the frequency I ran into a strange bug while implementing a complex band-pass filter. The filter design process is based on the windowing method and seems to be working well, but when i apply the coefs of the filter on my time domain input samples and perform DFFT the carrier frequency is doubled. I can´t really figure out why is this happening but my guess is that i'm doing something bad during the convolution. The first plot shows the magnitude response of the complex band pass filter. The second plot shows FFT of my input 5kHz carrier signal, the third one shows the same data but with BP filter applied. Everything is happening in the complex plain.

You complex multiplication code inputs 2 samples for every output sample.  That will act as a frequency doubler.


Why is deconvolution used in telescopes, CT scans and seismographs? Can you explain to someone with no mathematical background why is deconvolution used in telescopes, CT scans and seismographs? How do you explain that the same technique can be used for different things. The distortions in telescopes, CT scans and seismographs are all different from each other.

As long as the "distortion" is caused by a linear system, its effect can be modeled by a convolution operation with some "filter". The filter is different for different applications and perhaps even different scenarios and is usually measured / estimated before the deconvolution operation.

The deconvolution itself attempts to recover the original signal given the measured output which is the convolution of a signal with a filter (essentially the signal was filtered and we use the result to try and recover the original signal before filtering).


Estimating Plane Pose Without Knowledge of Intrinsic Camera Parameters Assume a camera has no skew, and has square pixels, that is, its camera-calibration matrix, $K$, looks like:
$$
K = \left[\begin{matrix} \alpha&0&u_x\\0&\alpha&u_y\\0&0&1 \end{matrix}\right]
$$

Intuitively, it seems that it should be possible to recover the pose of a rectangular planar object, $P$, (for example a flat piece of paper) by only knowing the projected image coordinates $\vec{p_i}$ for $i={1, 2, 3, 4}$ of the object's four corners without knowledge of $K$.  To be more precise, I only care about the orientation of the unit normal vector, $\vec{n}$, of $P$, not its location.

Following the answer given in response to this question:  Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers,

We choose world-coordinates such that $P$ is the plane $Z=0$.  We find the homography, $H$, mapping the four-corners of $P$ in world-frame coordinates to its projected coordinates $\vec{p_i}$ for $i=1,2,3,4$.

Then, $H = K[R | \vec{t}]$.  Where $R \in SO(3)$ is the camera's rotation and $t \in \mathbb{R}^3$ is the camera's translation.  If we make the assumption that $K$ is identity (I'll return to this later), then
$H = [\vec{R_1}|\vec{R_2}|\vec{t}]$.  Where $\vec{R_1}, \vec{R_2}$ are the first and second columns of $R$.  The third column of $R$,  $\vec{R_3}$, is equal to $\vec{R_1}\times\vec{R_2}$.  Finally, $\vec{n} = R\left[\begin{matrix}0\\0\\1\end{matrix}\right] = \vec{R_1}\times\vec{R_2}$ and we are done.

Now let's examine the effect of $K$ on our answer $\vec{n}$.  My intuition tells me that $K$ shouldn't matter because varying $\alpha$ will ultimately act as a scalar multiple on our answer $\vec{n}$ which is normalized out since I am only concerned with the direction of $\vec{n}$, not its magnitude.  In addition, it seems that varying $u_x, u_y$ should have the effect of translating the locations of the projected corners $\vec{p_i}$, which should affect $\vec{t}$, but not $\vec{n}$.

Let's see if this is true.  If $K \neq I$, then $H = K[\vec{R_1}|\vec{R_2}|\vec{t}]$ and $$\begin{equation}
\vec{n} = K \vec{R_1} \times K \vec{R_2} = det(K)K^{-T}(\vec{R_1}\times\vec{R_2})\end{equation}$$

In particular, if $\alpha \neq 1$,
$$\vec{n} = \alpha(\vec{R_1} \times \vec{R_2})$$
which is our answer from above times a scalar as expected.

On the other hand, if $u_x, u_y \neq 0$ then
$$det(K)K^{-T} = \left[\begin{matrix} 1&0&0\\0&1&0\\-u_x&-u_y&1 \end{matrix}\right]$$

and
$$\vec{n} = \left[\begin{matrix} 1&0&0\\0&1&0\\-u_x&-u_y&1 \end{matrix}\right](\vec{R_1} \times \vec{R_2})$$

In other words, our choice of camera center has a large effect on our calculations.  This is not intuitive to me.  Isn't our choice of camera center, $u_x, u_y$ completely arbitrary (for example $[0, 0]$ vs. $[w/2, h/2]$) or is there something I am missing?

I think there is a mistake in those calculations.

How can you assume that $K$ has only a scaling effect? The principal point is integrated in $K$, it's not just a matrix of focal lengths.
How come is this true?

$$det(K)K^{-T} = \left[\begin{matrix} 1&0&0\\0&1&0\\-u_x&-u_y&1 \end{matrix}\right]$$

You don't end up with such equal numbers in the diagonal if you compute $det(K)K^{-T}$.

To me, your question is a bit unclear and I guess some assumptions are fishy.


Encounter 0 when calculating log power spectrum To convert power spectrum to a log-scaled one, how to define log10(X(k)) if X(k)=0 for some k?

For sake of illustration, I brief my process as follows which is a convention one:

Calculate the magnitude spectrum X(k) of the time-domain signal x(n) by X(k) = abs(fft(x(n))).
Convert the magnitude spectrum or power spectrum to db by 20*log10(X(k)) or 10*log10(X(k)**2), respectively.

My problem arises when there is X(k)=0 when computing log10(X(k)) which is either not defined or -inf. How to deal with this?

Typically, one's data has some sort of noise floor (-96 dB, etc.).  So one common way to deal with the FFT bins that are zero (or tiny) is to replace anything below the noise floor with the noise floor level, since any value below that value is most likely not useful data.  Doing this before taking the log() function may provide some computational efficiency.


Difference Between monocular video and stereo video? What is the difference between monocular video and stereo video ,as i am studying about motion tracking i need to know the meaning of these two things.

Monocular video is taken using a single camera.  Stereo video is taken using two cameras side-by-side, producing two video streams.  That's pretty much it.

A stereo video pair allows depth measurement, so you can find out how far away things are.  However, this is not an easy task for complex scenes.


Lowpass Lanczos filter I am new to Lanczos filters and trying to figure out how they can be used to generate a zero-phase lowpass filter with good time-domain properties.
This example shows how to implement such a filter in Python but does not give any justification for the choice of window length, despite it determining the filter order, which I would think a relatively important parameter. Also, I am not sure about its phase response.
Does anyone have insights about the uses of Lanczos filters in this sort of context, their pros and cons, and how to choose the window size objectively?

This is a rather late response, but I just came over this question. So if it could of any help to future problems for anyone, I will post an answer:

Although your question is pretty comprehensive; In my experience with the Lanczos filter, it is a rather versitale filter especially when it comes to image processing. I have mainly used it for scaling images, where the filter is a part of other FIR processes. To use it I would design the Sinc function(which is defined by the Lanczos resampling function) to be dependent of the following parameters: Upscale factor (N) and filter order (a). When I defined the window length it is a simple function of the two parameters which is:

k = (a * 2) + 1

and then the window space results in:

x = np.linspace(-a, a + (N-1), k*N)

where k*N is the length. The Lanczos Sinc function is then multiplied with the interval values.
This was to give a practical example. In other words, you need to know what you want to do with the filter. It is highly dependent on what you want to filter, and you need to define one (or more) parameter, besides the filter order, to define the actual window length to take full use of the filter.

To give a quick overview of the Lanczos:

It is just one of many ways to realize the Sinc function in practical terms, and it is quite good when it comes to image sharpening and interpolating between resampled discrete data points.
It does have some limitations though. If the range of values of the interpolated signal is wider than the range spanned by the discrete sample values, it can cause ringing artifacts which ultimately leads to clipping artifacts. However the effect of this is quite low compared to many other realizations of the Sinc function.


How do I find a 50Hz sin/cos signal in a measurement among other frequencies and noise? I measure magnetic fields with high precision, in a gradiometer setup with fluxgate magnetometers, sampling at 2000Hz. I can have big data sets if I want to - one is 600000+ samples long. I need to find the phase and amplitude of a 50Hz sine signal. However there are other (low frequency - 2.145Hz) stronger signals superimposed, along with noise. Because of those low frequency signals the DC offset is hard to eliminate, too. As a result, an FFT gives no results yet.

Now it occured to me that it might be easier to find the 50Hz signal I look for (which is several orders of magnitude smaller then the DC offset and the other 2.1Hz signal) and determine its phase and amplitude, in stead of filtering out the rest.

What are good approaches and algorithms for that? I code this in python.

edit: here is my python (with numpy) code for my carlos loop algorithm.

   def decode_with_carlos_loop(self, frequency, periods_num):
        offsetfree_data = self.file_data
        shape = offsetfree_data.shape
        sampling_freq = 2000

        chunk_size = int(sampling_freq/frequency)
        alpha_step = .075
        indices = np.arange(chunk_size)

        for sensor_cnt in range(0, 4):
            for axis_cnt in range(1,3):
                error_sin = error_cos = 0
                cos_lo = sin_lo = []
                theta_last = 0
                theta = np.zeros(int(shape[1] / chunk_size) + 1)

                for chunk_cnt, chunk_pointer in enumerate(range(0,shape[1],chunk_size)):
                    chunk = offsetfree_data[sensor_cnt, chunk_pointer:chunk_pointer+chunk_size, axis_cnt]

                    cos_lo = np.cos(2*np.pi*frequency*indices/sampling_freq + theta[chunk_cnt])
                    sin_lo = np.sin(2*np.pi*frequency*indices/sampling_freq + theta[chunk_cnt])

                    # mean is the lowpass filter
                    error_sin = np.mean(np.multiply(chunk, sin_lo))
                    error_cos = np.mean(np.multiply(chunk, cos_lo))
                    error = error_cos*error_sin
                    theta_last = theta[chunk_cnt+1] = theta[chunk_cnt] - alpha_step*error
                print("sensor:axis:" + str(sensor_cnt) + ":" + str(axis_cnt) + "\tlast theta: " + str((theta_last + np.pi) % (2 * np.pi) - np.pi) + "\tcos_lo: " + str(cos_lo[0]) + "\tsin_lo: " + str(sin_lo[0]) )

You can try tracking Phase-frequency of your $50 Hz$ signal using Costas Loop. Costas loop does not require the signal to be pre-processed in order to expose the desired frequency.

I am not giving details of a Costas Loop because it can be found anywhere.It is pretty popular Carrier Recovery technique and a good starting point would be Wikipedia: CostasLoop

Instead, I would like to tell why I have chosen Costas Loop over other more common Squared-Diffrence Loop and Phase-Locked Loop :

Squared-Difference Loop and PLL requires the signal to be pre-processed by a squaring non-linearity and Band-Pass filtering at $f_{center} = 2*f_c$. This is done to emphasize the desired frequency component. This preprocessing step is not required for Costas Loop.
Phase Error sensitivity of Costas Loop is approximately double compared to PLL and Squared-Difference Loop. So, even smaller phase offsets in desired Frequencies are more accurately locked onto.
There is only Low-Pass filtering required in Costas Loop which can be implemented pretty easily as a Moving Average implementation. As DC-Offset and other low frequency noise is very high which is leaving FFT technique not useful, you can use a low cutoff LPF of enough taps to get a sharp transition, in order to get accurate and almost noise free Error signal $cos(2(\phi - \theta))$.

I have the MATLAB code implemented and customized for your situation at below path:

CostasLoopPhaseRecovery

The MATLAB Code runs n_runs number of Monty-carlo simulations to show that the algorithm will converge to the true phase of desired frequency eventually.

There are few design parameters which will depend on scenario to scenario. Like, in your case since there is a Large DC offset and a very low frequency component of large magnitude, therefore, you will have to use a good Low-Pass Filter to filter out the Phase Difference signal $cos(2(\phi - \theta))$. I have used a Moving Average Filter for the purpose of Low-Pass Filtering, and so I had to increase the length of filter and increase step size $\alpha$ in order to get phase convergence accurate and faster.

You will see a plot of Phase Convergence looking like below:

The MATLAB code assumes a DC Offset muchlarger in magnitude than the desired signal component amplitude. I have added $\phi = 0.2$ in the desired signal and the Costas Phase recovery Loop converges to $\phi = 0.2$. There is an inherent ambiguity in phase recovery of $\hat{\phi} = \theta + n\pi$, which appears in the plot too, and it depends on the initial phase of the desired frequency which locally generated and multiplied to the incoming signal.

Amplitude Estimation:

Once you have a pretty good estimate of the phase $\phi$ of the sinusoid at frequency $50Hz$ you can generate a reference signal $x[n] = cos(2\pi. 50.nT_s + \hat{\phi}), \forall n \in \{ 0,1,2,3, \cdots , N-1\}$, where A is the parameter to be estimated. You can now use Least Squares technique to estimate the Amplitude as follows:
$$\hat{A} = \frac{1}{N} <x,y>$$
where, $<x,y>$ denotes inner product. The problem is that in doing so we have ignored the fact that our Noise is not uncorrelated (or White) but colored. So, this might result into very wrong estimates. The way to fix the estimation is by generating reference for whatever known frequencies are there in your signal and estimating Amplitude of them, and modeling other unknown frequencies as Colored Noise.

So, in your case, you know that there is large DC offset in your signal and some small frequency component around $2.5Hz$. Assume DC, $2.5Hz$Sinusoid and $50Hz$Sinusoid amplitude as $A_o, A_{2.5} and A_{50}$. Let $y[n]$ be the measured signal, and then you can model $y[n]$ signal as : $$A_o.cos(2\pi 0.nT_s) + A_{2.5}cos(2\pi 2.5nT_s) + A_{50}cos(2\pi 50 nT_s + \hat{\phi}) + w(nT_s),$$where $w(nTs)$ is sampled colored noise (meaning Correlated).

In Matrix Form it would be:
$$\begin{pmatrix} y \end{pmatrix} = \begin{pmatrix} cos(2\pi 0.0T_s) & cos(2\pi 2.5.0T_s) & cos(2\pi 50.0T_s + \hat{\phi})\\ cos(2\pi 0.1T_s) & cos(2\pi 2.5.1T_s) & cos(2\pi 50.1T_s + \hat{\phi}) \\ \vdots&\vdots&\vdots\\cos(2\pi 0.N-1T_s) & cos(2\pi 2.5.N-1T_s) & cos(2\pi 50.N-1T_s + \hat{\phi})\end{pmatrix}. \begin{pmatrix}A_0\\A_{2.5}\\A_{50} \end{pmatrix} + \begin{pmatrix} w \end{pmatrix}$$ $$y = S.A+w$$
Then the LS solution would be the following:

$$\hat{A} = (S^HS)^{-1}S^H.y$$You can see that the noise is ignored again even though it is colored. The way to fix this is by estimating the Noise Covariance matrix and decolorising the noise or Whitening the noise, and then applying LS Technique.

You can also read about MAFI Algorithm to estimate Amplitudes of known Sinusoids in Colored Noise. That I hope will help you definitely. MAFI runs pretty close to Cramer-Rao Bound even at Low SNRs.


regarding visual dsp ++ i've written a source code in visual dp ++ in c language. the platform is sharc processor adsp 21062.after the build is complete. while running the program it is displaying a message
 Instruction timed out with PC at: 0x20a5a . what may be the error in my c code. can you please resolve?

The processor is stalled. For example you are trying to read a peripheral register and there is no data available. This can have many different reasons so without further detail that's impossible to assess. First step would be to look at program memory at address 0x20a5a and see  what the processor was trying to do. Take a look at the disassembly window.


What is the technical term for the product of "number of pixels" and "number of channel per pixel"? Assume I have m x n pixels of an RGB image. If I calculate 3mn, what is the name of quantity for this value?

number of channel per pixel

This is typically referred to as colour depth or bit depth and expressed as bits or 'bits per pixel' (bpp)

  number of pixels

There is no special terminology for this, it's just "Pixel count" or similar. It's unusual to need to refer to the total pixel count of an image without any other context, the dimensions are usually used instead.

  If I calculate 3mn, what is the name of quantity for this value?

There is no special word for this quantity, but the units would typically be bits or bytes, for example

1024 x 768 x 32bpp

= 3,145,728 bytes (3Mb)


Expressing 2N point DFT in terms of N point DFT I have a problem with expressing odd samples of X2 in terms of X1. I understand that the resulting DFT will be more precise in terms of expressing the exact spectrum of signal x[n], due to more samples. Moreover I know that even samples of X2(k) are copy of the spectrum of X1(k), however I do not know how to mathematically compute the odd ones.

The task is as follows:

Let $x[n]$  be  a  periodic  sequence  with  period $N_1$ .   Thus $x[n]$  is  also  periodic  for  period $N_2 =  2 N_1$ .   We  may compute $X_1[k]$:  $N$-point DFT of $x[n]$ and $X_2[k]$:  $2N$-point DFT of $x[n]$.

Express $X_2$in terms of $X_1$ Hint: it is easy with even samples $X_2[2m]$ , harder with odd ones $X2_[2m+ 1]$ , for $m$ integer.

Consider a sequence $x[n]$ of length $N$ whose $N$-point DFT is $X[k]$. Then let $X_2[k]$ be the $2N$-point DFT of $x[n]$.

As you have stated, the even indexed samples of $X_2[k]$ will be easily shown to be: $$X_2[k] = X[k/2] ~~~,~~~k = 2m, m=0,1,...,N-1$$

Then the odd indexed samples of $X_2[k]$ will be given by the $N$-point DFT of the signal  $x[n] e^{-j \frac{\pi}{N} n } $, $n=0,1,...,N-1$.

$$X_{2}[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{\pi}{N} n} e^{-j \frac{2\pi}{N} k n}  = X[k + 0.5]  ~~~~ ,~~~ k =2m+1, m=0,1,...,N-1 $$


In Convolutional Neural Nets, what do convolutions look like? The early stages of the Convolutional Neural Networks are performing classical convolutions with a certain kernel size on the input image. Is is possible to express in common terms the type of convolutions they are generally performing (lowpass, highpass, differential, anisotropic...), or are they completely application-dependent/impossible to interpret ?

Alternatively, can you show some sample kernel ? I would be very grateful.

The kernels used by a ConvNet are nothing but neural weights. You can think of them as a multilayer perceptron with some connections cut off and some weights restricted to be equal (weight sharing).

With this in mind, we must take into account that the kernels (or filters, in this context) are learned, so they depend exclusively on the type of inputs and the desired outputs. Each learning process is unique for a set of inputs/outputs, and so there is no way to know a priori what the filters will turn out to be. They are, as you state in your question, application-dependant.

They can be interpreted, though. Take for example a simple case, where we are trying to classify images that consist of horizontal lines only, and others of vertical lines only. If you do this, then the kernels will indeed look like that: vertical and horizontal lines. That's because the filters have learned common patterns to look for on the inputs, and they have acquired that form in order to be able to match those patterns and maximize their correlation to them.

In Figure 3 at Krizhevsky et al., there are some kernels for you to see. Kernels looking for horizontal lines will consist of horizontal lines. Kernels looking for vertical lines will consist of vertical lines. And so on.

To see this more clearly, you can train a CNN and take one of the kernels learned and filter any image with it to see what happens.

I think that a paragraph from this website can be of relevance:

  So what can we conclude from these feature maps? It's clear there is spatial structure here beyond what we'd expect at random: many of the features have clear sub-regions of light and dark. That shows our network really is learning things related to the spatial structure. However, beyond that, it's difficult to see what these feature detectors are learning. Certainly, we're not learning (say) the Gabor filters which have been used in many traditional approaches to image recognition. In fact, there's now a lot of work on better understanding the features learnt by convolutional networks. If you're interested in following up on that work, I suggest starting with the paper Visualizing and Understanding Convolutional Networks by Matthew Zeiler and Rob Fergus (2013).

In the paper mentioned in the previous quote, there are lots of examples of learned kernels, related to the images that were used to train the network. I think that paper can give you some useful insight.


Nutrition facts label OCR My team is working on a program that can read nutrition facts labels using OCR, much like NutriScanner (though by its reviews it doesn't seem like it works particularly well).

I've seen other questions noting that it'd be a good idea to straighten out the image prior to running OCR on it. Is there a good image library that can automatically do that? I've looked around and found numerous Photoshop tutorials on how to do it manually, but that's not what I'm looking for.

I know OCR engines like ABBYY have some pre-processing features built-in, but I'd prefer to piece together a solution using Tesseract and a free library for the pre-processing.

I haven't found much in the way of leveraging the positioning of the items on the label to improve scanning accuracy, but any suggestions would be appreciated.

Been there, done that.  One of our on-going projects is to create a process for processing nutrition labels from iPhone camera for a health-tracking app.  My company decided to develop a solution for this particular application - extracting data from US nutrition labels.  This solution will be used for this client, but we decided to pre-package a few other flexible capabilities along the way, to be used by a wider audience.  The solution will be available to general public in about 3 weeks.  (I'll come back here and post an update.)

Main goal that will make or break the entire idea is to get usable images.  Curved images, heavy shading, low resolution, and blur will all substantially decrease OCR quality, often to no quality.  In apps, we found that user training and guidance is the most helpful method.  Then there are some technical tools such as quickly detecting quality of image and suggesting to re-take it if needed.  In general, if you can achieve "high quality fuel" (images) you can expect high performance from your "machine".

Next, the OCR.  Tesseract can do good job reading clear and simple test such as pages from books.  Commercial products such as ABBYY distinguish themselves by working well on tough images - shading, distortions, small prints, etc.  And unfortunately if mobile images are used, they are more often on bad side than good quality side.

Next, pick your approach to locating and extracting data.  Please see my answer on text parsing vs specialized test extraction tools:
https://stackoverflow.com/questions/3070732/processing-ocred-text
For our project we'll be using ABBYY FlexiCapture for targeted nutrition data.  It has special tools even when OCR makes mistakes to still find appropriate data (sort of controlled fuzzy search), and that was an important factor in selecting it for the task.


An alternative algorithm for blob detection Here's a sample image of my experiment which is to be processed for finding the information of all the blobs (centre, dimension etc.). Basically the image contains some noise, few bold lines, and of course few blobs.

Full size image (12Mb)

My aim is to program a FPGA using an image processing package "Sapera APF". Sapera APF is an application for programming Teledyne DALSA frame grabbers to perform image processing. And unfortunately it doesn't allow us to write code in conventional way, neither allows to use basic coding features (if else statements, loops etc.) but only supports some inbuilt functions in form of code blocks, and we can only arrange these blocks.
FYI, Here's a link to a list of functions supported by Sapera:
https://www.dropbox.com/s/9f5qc9uki0lsalu/sapera_functions.docx

Now my first objective was to test my algorithm on matlab thereafter implementing it on Sapera. Here is my approach in matlab:

Applied mean filter (for step 3)
Thresholded image using adaptive thresholding (returning a binary image)
Removed lines by checking connectivity, all lines are connected with above 100 pixels      (for better connectivity mean filter was applied first)
Now since only blobs were remaining, their locations were identified (with grayvalue=255) and corresponding pixels were taken in original image.

This is a simple algorithm which works fine in matlab. But the real problem is that for step 2 (i.e. connectivity) of course I had to write loops going to each pixel checking for connectivity, which can't be done in Sapera. So I am looking for some other simpler algorithm to separate blobs in my images, without using loops and conditional statement, so that it can be implemented in Sapera, but now I just need its matlab code or algorithm.

Edit:
I cannot crop out strip containing blobs because other images may contain lines very close (maybe overlapping) to blob. For example, https://www.dropbox.com/s/m2j5pqo6stg4s1z/h125_s_-30.bmp

Since you have convolution and dilation you can use the Laplacian of Gaussian blob detector to detect the blobs and dilation to find the centre points. In MATLAB:

sigma = 3;

% Filter image with LoG
I = double(rgb2gray(imread('NCxQ8.jpg')));
h = fspecial('log', sigma*6, sigma);
B = -imfilter(I,h);
imagesc(B); pause;

% Threshold
B(B < 5) = 0;
imagesc(B); pause;

% Dilate then compare to get all
% local maximum in 3 pixel radius
Bd = imdilate(B,strel('disk',3));
imagesc(Bd); pause;
P = (B == Bd) .* B;
imagesc(P);

% Get point locations
loc = find(P > 0);
fprintf('There are %d blobs\n',numel(loc));

Note: "NCxQ8.jpg" is the smaller image embedded in the post. Adjust sigma for larger ones.


Are there any efficient algorithm to detect circular shape of arbitrary size in an image? I have images that are 1278 by 1024. They are pictures of a certain product that is rectangular in shape and are about 25mm by 20mm. As a part of a quality inspection process using computer vision, the goal is to detect contamination. The contamination are roughly circular in shape (caused by splattering of a certain material). Due to the special lighting, the contamination is seen as a roughly circular white area, surrounded by a green area (also circular, ring-like, concentric to the white area).

So basically it looks like a green donut, with a white 'hole'.

The goal is to detect such defect. However, the size of the contamination can vary. Furthermore, the size of the 'hole' can vary even if the outer diameter is the same.

I think I need to define a minimum and maximum size that I look for, and scan the image, but more advice and detail is needed.

Thanks

Usually this problem is tackled using the Hough transform. Check it here:

OpenCV
Scikit

After you have some hints for the possible location of the circles and their diameter then you can calculate the distance (e.g. Euclidean) from the pixels you have to an ideal model (circle parametric equations: x = r.cos a; y = r.cos a)


Is sparsity induced penalty in source separation "Entrywise matrix norms"? I am reading this paper where they introduce norm penalties for source separation. In table 1, the $\log/ l_1$ type is $\sum_{g} log(\epsilon +  \lVert H_{g} \rVert_1)$. I wonder this $\lVert H_{g} \rVert_1$ refer to which definition of $l1$ norm? In wikipedia, it has 2 definitions for this symbol.

"Entrywise" matrix norms: $\Vert A \Vert_p = \Vert \mathrm{vec}(A) \Vert_p = \left( \sum_{i=1}^m \sum_{j=1}^n |a_{ij}|^p \right)^{1/p}$ (wiki) or
Matrix norms induced by vector norms: $ \|A\|_1 = \max_{1 \leq j \leq n} \sum_{i=1}^m | a_{ij} |$ (wiki)

which is simply the maximum absolute column sum of the matrix.

One thing is below the table, the description said ""all norms are elementwise over matrix entries". So I doubt it is the first definition but not fully understand. Very appreciate if you may add explanation why it is that type of norm.

I haven't read the paper, but from the information you provided it appears the appropriate definition is the first one you gave, i.e., the "entrywise" matrix norm.

The text you quoted refers to the norms being "elementwise" which is equivalent to being "entrywise". One way to look at this norm is as the $l_1$ equivalent of the Frobenius norm. Essentially, you are replacing the sum of the squares of the matrix elements with the sum of the absolute value of the elements.

Trying to  minimize the $||H_g||_1 $ norm using this definition, is an approximation for trying to reduce the number of non-zero entries in the matrix $H_g$. You would have to examine the problem formulation in the paper to see if that makes sense.


Problem with scipy.signal.butter for designing filter bank I need to design a filter bank of band pass filters (code for implementation is shown below). There is a problem in that if my center frequencies fc are much smaller than the sampling frequency fs then you get the behavior as shown in the first graph:

A temporary solution is to change the sampling frequency to fs = fh*4 where 'fh' is the hi cutoff frequency of each respective band. This gives

Which is what I want, however because the sampling frequency is no longer the true one (relative to the time series I then want to filter) I cant use it for filtering purposes using lfilter(b,a,x) where x is my time series because all of the bands will be identical since they have the same non-dimensional frequency! My question is what can I do to make this work as expected for a true fs? Presumably there is an issue with finite precision?

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

def butter_bandpass(lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def fbank(fs):
    """
    Calculate third octave filter bank.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    # ISO standard center frequencies (bins)
    F0 = [50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500,
          630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000,
          5000, 6300, 8000, 10000, 12500, 16000, 20000]

    plt.figure(figsize=(10,6))
    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        fs = fh*4 # Temporary solution
        fs = 5e4
        b, a = butter_bandpass(fl,fh,fs)
        w, h = freqz(b, a, worN=2000)
        plt.loglog((fs * 0.5 / np.pi) * w, abs(h), label="fc = %d" % fc)

    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Gain")
    plt.legend(loc='best')
    plt.savefig("test_filter.png", bbox_inches='tight')
    plt.show()
    return 0

Although Im not sure how I can display the frequency response using this method, I have found that it is workable to use second-order section method for filtering. So say we have a signal x with sampling frequency fs_real that needs to filtered using a constant percentage bandwidth (in this case third octaves) then the corresponding power spectrum will be given by the get_cpb function:

def butter_bandpass_sos(lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = butter(order, [low, high], btype='band', output='sos')
    return sos

def get_cpb(x,fs_real):
    """
    Calculate third octave spectrum from pressure time series.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    cpb = []

    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        sos = butter_bandpass_sos(fl,fh,fs_real)
        y = sosfilt(sos,x)
        cpb.append( 20*np.log10(np.std(y)/2e-5) )

    return F0, cpb


Any simple Digital All-Pass Filter Design Tools? I would like to be able to implement an All-Pass filter with a frequency-dependent group delay.  I need a maximum group delay at low frequency of about 20 samples (for Fs = 44.1kHz), and it needs to fall to zero(-ish) at the Nyquist frequency.  Ideally, I would like to specify a corner frequency and a rate at which the delay falls off either side of it.  How the heck would I design such a thing?

I'm afraid I'm not sufficiently competent to follow the design approaches I can find on the web that plunge straight into poles, zeros, transfer functions and complex conjugates.  What I can do, though, is implement the darned thing efficiently in code once it is designed.

I am able to use Matlab to design FIR and IIR low-pass filters, because it has a built-in tool that spits them out.  But it doesn't have an All-Pass filter design tool (at least not one I can find).

This is an excerpt from the FilterScript user's guide, which should provide you with enough information on how to design the filter. Please refer to the other product resources on the tool's homepage for more information.


STFT computational burden I have a signal with length of $N$.
The sampling frequency is $f_s$.
Window size is $w$
The noverlap size is set to $m$.
And $nfft$ is selected for frequency resolutions.

I am running a STFT with matlab spectrogram function.

[S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)

what is the effect of nfft , $f_s$ and  noverlap on STFT
computational burden?
what is computation complexity of STFT with these parameters size?
{in term of big O() }

Note: The FFT computational complexity for a signal with length of $N$ is $Nlog(N)$

The number of FFTs is roughly N/(w-m). An fft costs you nfft*log2(nfft) so the total is proportional to N/(w-m)*nfft*log2(nfft).

A typical application would make nfft = w and m = w/2, so this would come out to be proportional to N*log2(nfft). So it's proportional to the length or you data set (obviously) and the log of the FFT size.


White noise and entropy confusion about concepts The spectral entropy (SE) of a signal is a measure of its spectral power distribution. The concept is based on the Shannon entropy, or information entropy, in information theory. The SE treats the signal's normalized power distribution in the frequency domain as a probability distribution, and calculates the Shannon entropy of it.

In general, Shannons' entropy of a signal is in between 0 and 1. 0 means no information content is present in the signal and 1 means otherwise. In communications, if I recall correctly we are uncertain of what the next symbol should be and so we prefer uncertainty as it conveys more information. My question is what if the entropy of the signal (spectral entropy) is say 0.12 i.e., quite small. For example, in a communication channel there is white noise but as soon as somebody speaks, the spectral entropy decreases (please see the example of detecting sine wave from white noise https://www.mathworks.com/help/signal/ref/pentropy.html)

Therefore, the information content of the signal now reduces during the time interval of the duration of the speech signal.

QUESTION: Does the reduced information content of the signal imply that the speech signal contains no information? How come a speech signal has no information or worth? Do we prefer high or low entropy (spectral and Shannon entropy)?

My confusions are:

High entropy --> more information content -- is that preferred or is useless?

Low entropy --> less information content -- Is that preferred?

Please correct me where wrong.

A signal's spectral entropy is not really a measure of information. Rather, it is a statistical measure of dispersion in the signal's power spectrum.

In other words, a random noise signal will tend to have large spectral entropy because its power spectrum is distributed among all frequencies. Speech signals tend to concentrate their power in a small number of bins, and their power spectrums are not flat. This results in lower spectral entropy.

A few comments about entropy:

Signals have information; sources have entropy. Signals never have entropy, and any reference to such concept is an abuse of terminology.
In communications, we prefer to have sources with large entropy, because then the messages that we transmit will convey a large amount of information. If you have a source with low entropy, then you should apply a source encoder to turn it into a large entropy source.
Shannon's entropy deals with uncorrelated sources. Many real-world sources produce correlated signals (such as text or speech) and their entropy is harder to define and calculate.


What's wrong with my Goertzel algorithm implementation? I've implemented Goertzel algorithm, now very simply in Python. But I cannot obtain the correct answer it's supposed to produce for a single DTF frequency. here is my code in Python:

import numpy as np

# omega frequency: radians/sample
omega = 0.23

# length of the data to generate
length = 200

# updated parameters
dft = 0.0
s1 = 0.0
s2 = 0.0

# initial phase
radians = 0.0
for i in range(length):

    # generate signal value
    value = np.cos(radians)

    # goertzel algorithm update
    s0 = value + 2.0 * np.cos(omega) * s1  - s2
    s2 = s1
    s1 = s0

    # dft sum update
    dft += value * np.exp(-1j * omega * i)

    # next value
    radians += omega

# final goertzel algorithm update
s0 = 2.0 * s1 * np.cos(omega) - s2
goe = s0 - np.exp(-1j * omega) * s1

# print results
print "DFT: {:.4f}".format(dft / length)
print "GOE: {:.4f}".format(goe / length)

Output:

DFT: 0.4979-0.0097j
GOE: -0.2065+0.4531j

The correct answer (with infinite data) is 0.5 + 0.0j.

DFT is therefore providing the correct result.

What's wrong with my Goertzel algorithm implementation?

https://en.wikipedia.org/wiki/Goertzel_algorithm

Why do you think the two should be equal?

Check out the final line of equation (6) in the article you link to:

I believe your misunderstanding is to not take account of the $e^{j\omega_0n}$ term.


Spectral Entropy Calculation in MATLAB How do I calculate the Spectral Entropy of a signal in MATLAB ? I know the basic steps but it would be nice if someone can help,

Calculate the power spectrum of the signal using FFT command in MATLAB.
Calculate the Power Spectral Density using the power spectrum or using any other technique.
Normalize the Power Spectral Density  between $[0, 1]$, so that it can be treated as a probability density function $p_i $.
Calculate the Entropy $H(s) = -\sum p_i\log_2\left(p_i\right)$

Technically this is not a MATLAB-esque forum, but I can explain the steps in more detail for you: Suppose your input signal is $x[n]$, and its DFT is $X(f)$. For real signals you may use the one-sided DFT, since the other half would be redundant when you look at its Power Spectral Density. (PSD).

Once you compute the DFT of your signal, the PSD is simply $|X(f)|^2$. That is, you need to take the absolute magnitude of your DFT result, squared.

You now need to normalize the PSD such that it can be viewed as a Probability Density Function, (PDF). Thus, a normalized PSD, (let us call it $PSD_n$) will simply be:

$$
PSD_n(f) = \frac{PSD(f)}{\sum_{f=\frac{-fs}{2}}^{f = \frac{f_s}{2}} PSD(f)}
$$

Finally, your spectral entropy will be:

$$
E = -\sum_{f=\frac{-fs}{2}}^{f = \frac{f_s}{2}} PSD_n(f) log_2[PSD_n(f)]
$$


PSD Magnitude with Welch's Method I'm having some trouble getting consistent results with PSD calculations. If I just take a basic sine wave

fs = 20000;
len = 16384;
fsig = 1000;
Sig = 1.5 * sin( (1:len) * 2 * pi * fsig / fs);
fsig = 700;
Sig = Sig + 0.9 * sin( (1:len) * 2 * pi * fsig / fs);

and I run it through the follow two Welch's functions with MATLAB's pwelch:

[matWelchD1, matWelchF1] = pwelch(Sig, length(Sig), 0, [], fs);
[matWelchD2, matWelchF2] = pwelch(Sig, ceil(length(Sig/2)), 0, [], fs);

I get the $1\textrm{ kHz}$ peak at approximately $0.65$ and $0.26$. I would love to blame frequency leakage due to how the frequencies line up, and there is some support for this because the $0.65$ peak has a much narrower base, but I see a similar result on real world data with FFTs that have broadband frequencies.

Anyone know why the two implementations would be different?

Your value for $f_s$ is missing. Assuming you work above Nyquist sampling rate (i.e. $f_s > 2*f_{sig}$), I obtain the following results:

fs = 20000;
len = 16384;
fsig = 1000;
Sig = 1.5 * sin( (1:len) * 2 * pi * fsig / fs);
fsig = 700;
Sig = Sig + 0.9 * sin( (1:len) * 2 * pi * fsig / fs);

[matWelchD1, matWelchF1] = pwelch(Sig, length(Sig), 0, [], fs);
[matWelchD2, matWelchF2] = pwelch(Sig, ceil(length(Sig)/2), 0, [], fs);

subplot(2,1,1);
plot(matWelchF1, (matWelchD1), '-o'); grid;

subplot(2,1,2);
plot(matWelchF2, (matWelchD2), '-o'); grid;

df1 = matWelchF1(2) - matWelchF1(1);
df2 = matWelchF2(2) - matWelchF2(1);

df1
df2

Output:

df1 =

    1.2207

df2 =

    2.4414

As you can see, both frequencies are nicely represented in the figures.

One more comments on the window length (i.e. second parameter to pwelch): It determines how big the sections are, into which the signal is divided. Each section is multiplied with a Hamming window and then the FFT is taken. Afterwards, the values of all FFTs are summed together, yielding the PSD estimate.

I.e. putting the signal length as the window length, would result in a single section. For your stationary signal this is fine, but in reality you would want to have a smaller value, or leave the parameter out (or use the empty vector []), such that the PSD can be more accurately estimated.

As you can see, different window length create a different amount of frequency samples: The longer the window, the more frequency samples. NOte that pwelch calculates the Power Spectral Density, i.e. the power per Hertz. In order to get the energy of a frequency you need to multiply with the bandwidth of each bin:

>> sum(matWelchD1) * (matWelchF1(2)-matWelchF1(1))

ans =

    1.5300

>> sum(matWelchD2) * (matWelchF2(2)-matWelchF2(1))

ans =

    1.5300

Both PSDs contain the same overall energy. Since the frequency samples in the first PSD are closer to each other, the PSD has higher peaks for each frequency (to deliver the same overall power).


Correlation and Convolution Both Correlation and Convolution are displacement function, i.e they are used to slide the filter mask across the image.

Convolution is same as correlation except that the filter mask is rotated 180 degree before computing the sum of products.

We mostly use convolution. why it is widely used? why convolution is preferred over correlation?

We have,

 Convolution in time domain   = Multiplication in Frequency domain

Do we have any similar relation for correlation ?

Correlation in time domain  =  _______________ in Frequency domain ?

Convolution as associative, while correlation is not.  The associative property is what allows you to convolve multiple filters into a single one, and then convolve that one filter with the image.  That is the big advantage of convolution over correlation.


Bandpass filtering IQ samples from USRP I have an USRP (B205mini) with an RX frequency of 900 MHz and a bandwidth of 100 MHz. I want to distinguish between two raw, unmodulated signals at respectively 875 and 925 MHz.

My approach would be to bandpass filter the output of the USRP, but I am not quite sure how to design and implement filters based on baseband IQ samples. Can anyone provide references or explanations ?

There are many possible methods to isolate those 2 signals.  You can heterodyne your IQ signal up or down by 25MHz by multiplying by a complex exponential of plus or minus that offset frequency, then baseband filter each IQ component separately with either a real FIR or real IIR low-pass filter.  Or you can multiply a FIR low-pass filter kernel by those same complex exponentials.  Or you can strip the poles and zeros of any canonical IIR bandpass from either the bottom or top half of the complex plane, and convert that to complex difference equation form.
Or you can use an FFT/IFFT fast convolution filter algorithm with an offset (non-complex conjugate mirrored) response.  Or just use FFTs (of the complex IQ input) of moving windows to compare the negative FFT result frequency bins against the positive half bins.  The latter is a common source for SDR waterfall visualizations.

etc.


Under what conditions is the convolution of an input signal with the system's impulse response periodic? I'm currently solving the following convolution problem from Oppenheimer's book:

In the solution, it was stated that "$x(t)$ periodic implies $y(t)$ is periodic"
So I wondered if it's always the case that if $x(t)$ is periodic then so is the output from the convolution, I tried doing research about this in hope of finding any supporting evidence to this but I couldn't, and "circular convolution" is all around the place.

I am somehow rephrasinge what @Hilmar already answered. In the continuous-time setting, convolution cannot invent frequencies. Indeed, there is an intimate relationship between Fourier and convolution, that is not always taught properly. Recently, Michael Bronstein was Deriving convolution from first principles

Have you ever wondered what is so special about convolution? In this
post, I derive the convolution from first principles and show that it
naturally emerges from translational symmetry.

The basic assumptions that a system is linear, and invariant by shift or in time (LTI) imply the concept of convolution. And convolution implies that a sine undergoing an LTI system will still be a sine, albeit with a time-shift (phase change) and amplitude modification (it can even vanish). In other symbols, through an LTI:
$$ \sin(\omega t)\to A\sin(\omega t+\phi)$$
$A$ can be zero, so frequencies can disappear. Yet if  $\sin(\omega t)$ is not present in the original data (no $\omega$), it cannot reappear. Frequency creation from nothing is a feature of nonlinear or time-variant systems.


Is the first derivative operation on a signal a causal system? Please help me sort this issue out.

Consider a system whose output $y(t)$ is the first derivative of the input signal $x(t)$.

We can write the first derivative of an input signal as follows:

$$y(t)=\frac{dx(t)}{dt}=\lim_{\Delta t\to 0}\frac{x(t+\Delta t)-x(t)}{\Delta t}$$

Thus, we can see that the system's output must remember an infinitesimal change in future input in order to calculate the derivative, thus system possesses memory. Perfectly fine till here.

The problem comes in checking if it's causal. By the definition of the derivative above, the system has to subtract an infinitesimal future value with the current value, to calculate the derivative. Thus by definition, this has to be a non-causal system, right?

However, most of the textbooks and websites mention the first derivative as a causal system.

Please let me know if I am making any mistakes. Thank you.

If the derivative exists at the given point, then it doesn't matter if you look (infinitesimally) into the future or into the past, you can do both, because both will give the same result:

$$x'(t)=\lim_{h\to 0}\frac{x(t+h)-x(t)}{h}=\lim_{h\to 0}\frac{x(t)-x(t-h)}{h}\tag{1}$$

So a differentiator can be (theoretically) implemented by a causal system.

Clearly, the system is not stable, because a bounded input can cause an unbounded output. E.g., a step at the input causes an unbounded output signal. So if by realizable we understand causal and stable, the ideal differentiator is not realizable. However, good approximations in a certain frequency range can be realized with small errors, in the continuous domain as well as in the discrete domain.


Problems calculating 8-point FFT of an 8-point sine wave by hand I need a small help for my assignment, which is to solve the 8-point FFT without using Matlab or C.

Preliminary information:

8-point sine wave: [0, 0.7071, 1, 0.7071, 0, -0.7071, -1, -0.7071]

8-point complex exponential: [1, 0.7071-j0.7071, -j, -0.7071-j0.7071]

So here's the attempts I did:

1.) Bit reversing the 8-point sine wave and rearrange it to:

[0, 0, 1, -1, 0.7071, -0.7071, 0.7071, -0.7071]

2.) Perform each butterfly calculation:

First pass:
[0, 0, 0, 2, 0, 1.414, 0, 1.414]

Second pass:
[0, -j2, 0, 2, 0, 1.414-j1.414, 0, 1.414+j1.414]

And the last:
[0, -j4, 0, 2-j2, 0, 2.828, 0, 2.828-j2.828]

However, the last pass which is the result does not match with the actual FFT result in Matlab:

[0, -j4, 0, 0, 0, 0, 0, j4]

Any steps I might have missed? Opinions and pointers are greatly appreciated. Thanks.

I will start from very beginning. So having signal: $ x[n] = [0; \ 0.7071; \ 1; \ 0.7071; \ 0; \ -0.7071; \ -1; \ -0.7071] $ and given exponent: $ w[n] = [1; \ 0.7071 - 0.7071i; \ -i; \ -0.7071 - 0.7071i;\ -1;\ -0.7071 + 0.7071i; \ i; \ 0.7071 + 0.7071i]$.

First we need to inverse the bits and rearrange our signal. For 8 samples we will obtain following change in indices: $ [0;\ 1;\ 2;\ 3;\ 4;\ 5;\ 6;\ 7] \Rightarrow [0;\ 4;\ 2;\ 6;\ 1;\ 5;\ 3;\ 7] $, thus our signal becomes: $ x[n] = [0; \ 0; \ 1; \ -1; \ 0.7071; \ -0.7071; \ 0.7071; \ -0.7071] $.

Now according to following schematic:

We calculate "butterflies" in following manner:

Where: $ W_{N}^{k} = e^{\dfrac{-2\pi i k}{N}} $ (rotating vector on our complex plane with step angle of $2 \pi/N $ )

Starting from a first butterfly in first stage ($x[0]$ and $x[4]$ as the inputs) we have following:

$ A = 0, B = 0, W_8^0=w[0]=1$, giving output: $ [0; 0] $, next butterfly:

$ A = 1, B = -1, W_8^0=w[0]=1$, giving output: $ [0; 2] $, next butterfly:

$ A = 0.7071, B = -0.7071, W_8^0=w[0]=1$, giving output: $ [0; 1.4142] $, next butterfly is the same with output:  $ [0; 1.4142] $.

So now, after first pass our vector is as follows (same to the one obtained by Peter Griffin):

$x'[n]=[0; \ 0; \ 0; \ 2; \ 0; \ 1.4142; \ 0; \ 1.4142] $

At second stage, first butterfly consists of samples: $x'[0]$ and $x'[2] $, where $ W_8^0=1 $. We can therefore calculate the output as: $[0+1\cdot 0; \ 0-1\cdot 0] = [0; \ 0] $. These are outputs $x''[0]$ and $x''[2]$ of the second stage.

Next butterfly is given by inputs $x'[1]$ and $x'[3] $, where $ W_8^2=w[2]=-i $. We calculate output $x''[1]$, $x''[3]$ as: $[0-i\cdot 2; \ 0-(-i)\cdot 2] = [-2i; \ 2i] $

Two following butterflies are calculated in the same manner:

$[0+1\cdot 0; \ 0-1\cdot 0] = [0; \ 0] \Rightarrow x''[4]$, $x''[6]$

$[1.4142-i\cdot 1.4142; \ 1.4142-(-i)\cdot 1.4142] = [1.4142-1.4142i; \ 1.4142+1.4142 i] \Rightarrow x''[4], \ x''[6]$

So our vector after second stage is following:

$x''[n] = [0; \ -2i; \ 0; \ 2i; \ 0; \ 1.4142-1.4142i; \ 0; \ 1.4142+1.4142i] $ - you've made a mistake here.

Finally, we calculate results for third stage. By looking at the main schematic, first butterfly consists of samples: $x''[0], x''[4] = [0; 0]$. The $W_N^k$ coefficient is given by $W_8^0=w[0]=1$.

Calculating butterfly we obtain:

$[0+1\cdot 0; \ 0+(-1)\cdot 0] = [0; \ 0] $ - these are final FFT values $X[0], X[4]$.

Moving to next butterfly based on samples $x''[1], x''[5] = [-2i; 1.4142-1.4142i]$ with coefficient $W_8^1=w[1]=0.7071 - 0.7071i$. Result: $[-2i+(0.7071 - 0.7071i)\cdot (1.4142-1.4142i); \ -2i-(0.7071 - 0.7071i)\cdot (1.4142-1.4142i)] = [-4i; \ 0] $

And next butterfly based on samples $x''[2], x''[6] = [0; \ 0]$, with  $W_8^2=w[2]=-i$ will give $[0+(-i)0; \ 0-(-i)0] = [0; \ 0] $

And the last butter fly based on samples $x''[3], x''[7] = [2i; \ 1.4142+1.4142i]$, with $W_8^3=w[2]=-0.7071 - 0.7071i$ will result in: $[2i+(-0.7071 - 0.7071i)\cdot (1.4142+1.4142i); \ 2i-(-0.7071 - 0.7071i)\cdot (1.4142+1.4142i)] = [0; \ 4i] $

So our final vector is following:

$X[k] = [0; \ -4i; \ 0; \ 0; \ 0; \ 0; \ 0; \ 4i] $

Which means MATLAB is calculating FFT correctly ;)

One final remark - always take care about $W_N^k$ factor - look on the schematic and it's pretty obvious.

Good luck!


Can cross-c.(x,y) be positive if signals x and y do not share frequencies? The question is if the cross-correlation of two signals $x$ and $y$ can be non-zero even if $x$ and $y$ do not share frequencies.

My approach is the following:
The cross-correlation  $\rho_{xy}\left(\tau\right):=x\left(\tau\right)\star y\left(\tau\right)$  can be expressed in terms of a convolution:

$$
\begin{array}{rcl}
\rho_{xy}\left(\tau\right) & = & x\left(\tau\right)\star y\left(\tau\right)\\
 & = & \displaystyle\sum_{t=-\infty}^{^{\infty}}x\left(t\right)\bar{y}\left(t+\tau\right))\\
 & = & \left(\bar{x}\right)\left(-\tau\right)*y\left(\tau\right)
\end{array}$$

i.e. the cross-correlation of $ x\left(\tau\right) $ and $ y\left(\tau\right) $ equals the convolution of $ \bar{x}\left[-\tau\right] $ and $ y\left(\tau\right).$

The convolution theorem  $\mathcal{F}\left(x*y\right)=\mathcal{F}\left(x\right)\mathcal{F}\left(y\right) $ allows us to express the cross-correlation in terms of it's frequency content and replace the convolution by multiplication: $$ \mathcal{F}\left(\rho_{xy}\right)\left[\omega\right]=\mathcal{F}\left(\left(\bar{x}\right)\left(-\tau\right)*y\right) \left[\omega\right]=\bar{\mathcal{F}\left(x\right)}\left[\omega\right]\mathcal{F}\left(y\right)\left[\omega\right] $$
 Due to Parseval's Theorem, signal energy in frequency and time domain corresponds to: $$ \sum_{\omega=-\infty}^{\infty}\left|\mathcal{F}\left(\rho_{xy}\right)\left[\omega\right]\right|^{2}=\sum_{\tau=-\infty}^{\infty}\left|\rho_{xy}\left(\tau\right)\right|^{2}.$$

Hence the energy of the cross-correlation of the signals corresponds in time and frequency dimension. Thus, unless two different signals share frequencies, they cannot interfere. Two orthogonal signals, i.e. $\left\langle x,y\right\rangle =0 $, can have shared frequencies without having a non-zero cross-correlation.

The solution to the self-proposed problem was correct, as it stands.


Are the error probability of orthogonal signal and on/off signal the same? I have to write the error probability of orthogonal signal and on/off signal,and now i have finished the  error probability of on/off signal,but i don't know how to write the code about error probability of orthogonal signal.

Generally,we will introduce the orthogonal signal and on/off signal with the picture below ,first is orthogonal signal,second is on/off signal.

(a)  orthogonal signal

(b)  on/off signal

And the definition of orthogonal binary signal is that if we do the inner product with the signal 1 and signal 0 ,the result will become 0. I mean $\int s_0  s_1 dt = 0$

And i found that the on/off binary signal is also $\int s_0  s_1 dt = 0$,so can i say the on/off signal is the orthogonal binary signal too,so their error probability should be the same theoretically?

I found that the on/off binary signal is also $\int s_0  s_1 dt = 0$, so can I say the on/off signal is the orthogonal binary signal too, so their error probability should be the same theoretically?

No! That's a wrong inference: Just because the zero vector (the off signal in your case) is orthogonal to everything doesn't mean everything has the same error probability.

When comparing error probabilities, you first have to agree on a way of making the energy spent on each transmitted bit equal.
Otherwise, you might be comparing an On/off-keying transmitter with 1 kW to a BPSK transmitter with 1 mW output power – you'll notice that this comparison will probably be won by the on/off-keying, no matter whether the other one is better.

So, same average bit energy is what we demand. Since one symbol carries one bit in your example, we can set the average symbol energy to 1 for both cases (without loss of generality, as we can later scale), and compare then.

So, sadly, the book you're learning from doesn't do that energy normalization to begin with.

Let's say for the on-off case, $A_{OOK} = 1$, and $T_b=2$, so that the energy the on-period is $E_{On}=P_{On} T_b = A_{OOK}^2 T_b = 1^2 2 = 2$. Obviously, during the off-period, the power is $P_{Off}=0=E_{Off}$, so that the average energy is $E_{OOK}=1$ (you always assume you send as many 0s as 1s, otherwise your source coding is broken).

For what you call the "orthogonal case", the power is always the same – the sign of the signal doesn't matter to the power, since you square. So, we set $A_{ortho}=\frac1{\sqrt2}$ and get $P_{ortho}=A_{ortho}^2=\frac12$ and thus also $E_{ortho}=1$. So, now we're fair!

That means that the pictures you posted aren't fair: The $A$ from the "orthogonal signalling" picture should be $\sqrt2$ smaller than the one from OOK.

We also see that we need a "stronger" amplifier for the On/off case just to get the same signal power. That alone might be an argument why OOK isn't that good.

Now, the rest of the derivation of what an error probability is depends on whether you assume coherent or incoherent reception – you don't define that, but your textbook will certainly do.

You'll, in any case, notice that the two things have different symbol error rates at the same noise power when you normalize their power.


Signal bandwidth of an 802.15.4 signal at 2.4 GHz This signal uses Offset QPSK with a half sine pulse shaping filter for the purpose of  getting a constant signal envelop. Does anyone know the exact bandwidth of this signal after being pulse shaped, i.e. the nul to nul bandwidth of the main lobe?

I found this link that details the spectrum for 802.15.4: https://www.semanticscholar.org/paper/IEEE-802.15.4-PHY-analysis%3A-Power-spectrum-and-Gupta-Wilson/2892da39a7ecec1945595d46950438a5bb8777af

And pasted the plot below.

This is consistent with what I would expect for a half-sine pulse-shaped pulse for a symbol rate of 1 MHz (and that would apply to BPSK, QPSK, OQPSK etc where the power spectrum for a random data sequence is the Fourier Tranform squared of the base pulse that is transmitted.).  In Offset-QPSK the I and Q channels can be viewed as two independent (in data content) BPSK transmissions, in quadrature, offset in time by half a symbol. Thus each produces an independent but identical power spectrum that adds in power.

Below I compare the Fourier Transforms for a rectangular pulse of 1 us duration vs the half-sine windowed pulse of the same duration. I normalized the power of the two to be the same. Although the pulse shape increases the width of the main lobe, there is a clear spectral improvement by doing the half-sine pulse shape as demonstrated in this plot (Notably, the power spectrum of the rectangular pulse rolls of at a rate of $1/f^2$, or -20dB/decade, while the power spectrum of the half-sine pulse rolls of at a rate of $1/f^4$ in spectrum).

Fourier Transform of rectangular pulse (Red) vs half-sine pulse (Blue):

This clearly is not as good as that which can be achieved with more spectral efficient pulse shaping such as raised-cosine shapes but offers the significant advantage to low power systems of being a constant envelope modulation when used in an O-QPSK modulation. In this case we are able to drive the power amplifier (PA, which typically dominates power dissipation in the transceiver) into saturation without introducing AM/PM distortion. In saturation the PA can operate at a much higher efficiency (Tx power radiated at a given DC power in).


what is the starting point of a contour as given by cvConvexityDefect() function in OpenCV i have to distinguish between palm and elbow in my contour(I'm using cvconvexitydefect for that purpose). Here's the image:

I'm taking contour of left image.

the 2 topmost points in right image are shoulder(detected by other method)

bottom left 2 are darker than bottom right 2 (to distinguish between left and right hand)

using the code defect.start is darker than defect.end(but for left hand defect.start is palm while for right hand defect.start is elbow)

I have 2 questions:

So how does cvconvexity defect get it's defect.start?
In any case how to differentiate palm from elbow?

To be clear i thought that i could make out the palm from elbow just by thinking that start might be upper most point of a contour or leftmost point of that but it ain't so.

here's the relevant part of code:

cvFindContours( HeadLess1, mem, &contours, sizeof(CvContour), CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE, cvPoint(0,0));

for(ptr=contours;ptr!=0;ptr=ptr->h_next){
    tem=cvConvexHull2(ptr,0,1,0);
    tem=cvConvexityDefects(ptr,tem,mem);
    CvConvexityDefect defect,cd;
    defect.depth=0;

    for(int l=0;l<tem->total;l++){
        cd=*(CvConvexityDefect*)cvGetSeqElem(tem,l);
        if(cd.depth>defect.depth)defect=cd;
    }

    if((*(defect.start)).y>(*(defect.end)).y){
        elbow=(*(defect.start));palm=(*(defect.end));
        //a bad inference for the case when palm is above elbow
    } else {
        elbow=(*(defect.end));palm=(*(defect.start));
        //a bad inference same case
    }

    if(elbow.x<320)
    {
        q=100;
    }
    else
    {
        q=200;
    }

    cvCircle(HeadLess1,*(defect.start),10,cvScalarAll(q),CV_FILLED);
    cvCircle(HeadLess1, *(defect.end), 10, cvScalarAll(q+50), CV_FILLED);

}

Headless1 one is my image. Thanks in advance.

If you can obtain coordinates for the points you show in the right hand image, and you know the points relating to the shoulders, could you presume that the palms were further from the shoulder than the elbow is?

If you assume this to be true, then you simply need to find the distance from the shoulder to the elbow & palm and call the point with the larger distance the palm.

Of course, you could bend your arm in such a way that the palm would be closer to the elbow, so perhaps this isn't an ideal method...


Getting Vario (m/s) from a bmp180 I am building a variometer as a hobbie, this could be a duplicate from an existing question: How to use Kalman filter for altitude prediction based on barometer data? .

Details follow:

I have a https://es.wikipedia.org/wiki/NodeMCU running http://www.micropython.org/ , it works nice and I have attached a BMP180 barometric pressure sensor, so I can retrieve: temperature, pressure and altitude, however, while the sensor is fixed in my desk and if I read from three seconds I can get values that range from (982.409 to 983.956 ) in terms of altitude, see the log:  https://pastebin.com/C3RbuyDn . So there is an error margin of more than 1 meter .

I would like to use it as a variometer for paragliding, and I only have the pressure sensor, I read it is better to combile with a magnetic sensor, however I want to keep it simple, so it is possible to have an accurate variometer with just the pressure sensor and some algorithms behind and off course with no lag ( no more than 1 second ) ?

Is there any existent algorithm for getting the variation (m/s) with just one variable (altitute) and also considering this chip report the altitude ?

Please look at the data referenced in pastebin .

Thanks for your time.

If you add an accelerometer to the project, a Kalman filter can give a good estimation of vertical speed.

With only a barometric sensor, I don't think it's possible to reduce the lag below 1 second.

import numpy as np
import matplotlib.pyplot as plt
import random

from filterpy.kalman import KalmanFilter
from filterpy.common import Q_discrete_white_noise

random.seed(65537)

# Standard deviation of simulated sensor data
std_a = 0.075
std_h = 0.42

# Earth gravity
g_n = 9.80665

# Sensor sample period, s
dt = 0.05

kf1 = KalmanFilter(dim_x=2, dim_z=1)
kf1.H = np.array([[1, 0]])
kf1.F = np.array([[1, dt], [0, 1]])

kf2 = KalmanFilter(dim_x=3, dim_z=2)
kf2.H = np.array([[1, 0, 0], [0, 0, 1]])
kf2.F = np.array([[1, dt, dt * dt * 0.5], [0, 1, dt], [0, 0, 1]])

# initial process covariance
kf1.P = np.array([[1, 0], [0, 0.01]])
kf2.P = np.array([[std_h * std_h, 0, 0], [0, 0.0001, 0], [0, 0, std_a * std_a]])

# Process noise matrix
std = 0.0039
var = std * std
kf1.Q = Q_discrete_white_noise(dim=2, dt=dt, var=0.0625)
kf2.Q = Q_discrete_white_noise(dim=3, dt=dt, var=var)

# Measurement covariance
kf1.R *= np.array([[std_h * std_h]])
kf2.R *= np.array([[std_h * std_h, 0], [0, std_a * std_a]])

n = 300
r_n = 1.0 / n

t = np.zeros(n)
h_sim = np.zeros(n)
v_sim = np.zeros(n)
a_sim = np.zeros(n)
measured_h = np.zeros(n)
measured_a = np.zeros(n)
v_est_kf1 = np.zeros(n)
v_est_kf2 = np.zeros(n)

v0 = 0
v = v0
h = 0

for i in range(n):
    t[i] = (i * dt)
    a = 1.0 / 32 * g_n * np.sin(4 * np.pi * i * r_n)
    a_sim[i] = (a)
    v += a * dt
    v_sim[i] = v
    h += v * dt + (a * dt * dt) / 2
    h_sim[i] = h
    measured_a[i] = a + random.gauss(0, std_a)
    measured_h[i] = h + random.gauss(0, std_h)

# Compute the speed estimations

for i in range(n):
    kf1.predict()
    kf1.update(np.array([measured_h[i]]))
    v_est_kf1[i] = kf1.x[1]

    kf2.predict()
    kf2.update(np.array([[measured_h[i]], [measured_a[i]]]))
    v_est_kf2[i] = kf2.x[1]

# Plot the results

plt.figure(1, figsize=(8, 12), dpi=80)

plt.subplot(311)
plt.axis([0, n * dt, -0.75, 0.75])
plt.plot(t, measured_a, 'y+')
plt.plot(t, a_sim, 'r')
plt.title('Vertical acceleration - gravity, m/s^2')
plt.legend(('Measured', 'True (simulation)'), loc='best')

plt.subplot(312)
plt.axis([0, n * dt, -1, 7])
plt.plot(t, measured_h, 'c+')
plt.plot(t, h_sim, 'b')
plt.title('Altitude, m')
plt.legend(('Measured', 'True (simulation)'), loc='best')

plt.subplot(313)
plt.axis([0, n * dt, -0.5, 1.0])
plt.plot(t, v_sim, 'b')
plt.plot(t, v_est_kf1, 'r')
plt.plot(t, v_est_kf2, 'c')
plt.title('Vertical speed, m/s')
plt.legend(('True (simulation)', 'Kalman filter (1 sensor)', 'Kalman filter (2 sensors)'), loc='best')

plt.savefig('kalman_vsi.svg')
plt.show()

After some experiments with simulated sensor data, I wrote an Android app that uses smartphone sensors. Source code on GitHub


Impulse response of LTI Systems I am trying to solve the attached question. Please I would appreciate if some one can check my approach and tell me if I am correct and how can I proceed to find h0(t)

[![enter image description here][3]][3]

$h_0(t)$ is inverse Fourier transform of $H(jw)$. Your formula is OK, you can continue your calculation to practice your math manipulation, why not.
To check the result, you can remark that $H(jw)$ is rotated $sinc$ function. And $sinc$ must be the Fourier transform of a rectangular box function.


Queries on Kalman Filter I am trying to apply kalman filter for video processing , i am studying about it from different sources but it take me towards question that if i don't know that where my object come in frame mean i don't know about its angle (location in frame) then is it possible that Kalman filter automatically detect its measurements e.g location (as i already applied algorithms for specific object detection , so i only want that detected object tracking) and then track it in other frames ?

Other thing which come to my mind is that two objects (which i want to detect) come in the same frame and in the same measurements to is kalman filter take them as two different objects for tracking ?

Main thing is that in the first frame it will be required to select the object of detection because it is obvious that the algorithm will not know automatically which object you want to track if a number of things will be moving in the video scene.Lucas-Kannade method is one of the methods which can detect moving objects in a given video frame . If you do not specify object it will show movement of all moving objects.Thats why you have to tell the algorithm which thing you want to track. Second thing i want to tell you is about Kalman filter which i think you are interpreting wrong , what Kalman filter actually do is that it uses the system model and input signal and generates its own output which is actually a predicted value . Now using the observation from the system it calculate the correction to be applied to the predicted output and hence final results are better than observations obtained from the system. Hence very first step is that you need a system model and the measurements. In tracking algorithms generally, Image processing algorithm is applied to detect the object of interest and then Kalman filter is applied to make the detected location more precise. Do not think that it will predict on its own. I hope it will help you.


why does/ doesn't phase count when you analyze a voice signal? I am new to digital signal processing, I want to know why phase should not affect the analysis of a voice signal.

Phase will have an effect on the mathematical analysis of speech signals however the traditional practice of speech (audio) processing has evolved based on the observation that human listeners do not have a phase sensitive hearing mechanism, which indicates (supported by experiments) that phase may be neglected in most speech processing applications, at least when compared to the frequency magnitudes. You may search the site for other related answers. But I think the most extensive source will be the audio engineering society.


How the white and colored noise differ in time domain I know the white noise is uncorrelated. So a perfectly white signal has an autocorrelation equal to an impulse at zero. This implies that the corresponding signal in the time domain is any function with no correlation whatsoever between the different samples.

And I know that any colored noise can be thought of as filtered white noise, so this will still look like a somewhat random signal but with less "sharp" transitions between the samples because the high frequency components have been filtered out. This means that in colored noise, the low frequency noise will have higher variance than the higher frequency noise.

But other than that can you give more intuition of the difference between the white and colored noise in time domain?

Lets say I have two exact same signals and I add white noise to one and colored noise with the same total noise power to the other. Does that mean in the signal which I added colored noise the DC part of the signal experiences a noise with higher average variance than the AC part of the signal?

If I use a white noise of the same power instead of colored noise, then what do I lose and how bad this approximation would be in a system simulation?

White noise looks like the plot below on the left in the time domain. A key point is no matter how much you zoom into it, it will still generally look like this. Why? Because by definition every sample (let's assume a sampled system although this is not limited to that) is uncorrelated from the next in a white noise signal. The auto correlation is therefore an impulse.

If it happens to be a Gaussian White Noise process, then it will also have a magnitude distributed as shown in the histogram on the right, but these are two different characteristics, it can be White and not Gaussian, or Gaussian and not White.

Now for comparison, below is the same result for a "Colored Noise", and in this case, Gaussian distributed process.  Interestingly here the histogram of the distribution, given sufficient samples is identical, but the key point is as you zoom into the time domain plot on the left, we see the result of "memory" due to the bandwidth limiting that was done. Because of low pass filtering (in this case), the signal cannot transition from one sample to any other arbitrary location on the next sample; the location of the next sample depends on previous samples and is therefore constrained. (The more we low pass filter, the longer the effective memory depth). So there is more predictability on the next sample based on previous samples; samples are correlated to adjacent samples, and the auto-correlation function will no longer be an impulse but will be spread out with a width depending on how narrow the filtering is in frequency.

In the frequency domain, the power spectral density for the "Colored Noise" could look like the plot below; in particular this plot is the low pass filtered signal discussed above (but colored noise could be any selective band of frequencies, or even multiple bands). In contrast the white noise would have the same level of power spectral density for all frequencies. So in case this point wasn't obvious to all readers; white similar to white light that contains all frequencies vs colored with again the analogy to light in that individual colors are given by selective frequencies of light.