How to store the signal in a variable and plot the signal in MATLAB I am trying to store the signal I have created in a single variable and graph the signal in two subplots in one single graph, one for left channel, one for right channel.

Fs = 44100;
N = 400
t = (1:N)*(1/Fs);
freq = 500;

sound_array = zeros(N,2);
%starts N x 2 matrix
sound_array(:,1) = sin(2*pi*freq*t);
sound_array(:,2) = sin(2*pi*freq*t);

%Left and Right channels
sound_L = sound_array;
sound_L(:,2) = 0;
sound_R = sound_array;
sound_R(:,1) = 0;
plot(sound_array);

R = 400;
beep = [sound_array; zeros(R, 2)];

%Sound 1
sound(sound_L, Fs);
pause(1);
%Sound 2
sound(sound_R, Fs);
%Play Consecutive
sound(sound_array, Fs);

I suggest you to do something like (one plot with two curves):

figure
plot(sound_L(:,1), 'b'); % blue
hold on; grid on; % prevent from redrawing and add grid\
plot(sound_R(:,2), 'r'); % red
legend({'left', 'right'});

Or if you want two subplots:

figure
subplot(211);
plot(sound_L(:,1), 'b'); grid on;
subplot(212);
plot(sound_R(:,2), 'r'); grid on;

Please mind that in your case sound_L and sound_R are the same. You are not making any use of beep variable.


Word Error Rate weights Usually ASR systems are evaluated using WER (word error rate), which summarizes 3 types of changes when calculating the edit distance: insertions, deletions and substitutions. According to the wikipedia page, there are two versions:

A. each type is given the same score: 1, 1, 1 (S, D, I)

B. Hunt's version: 1, 0.5, 0.5

I've noticed that sclite has it's own weights

C. sclite: 4, 3, 3

What are the use cases for each weighting scheme? My goal is to compare different speech recognition APIs.

The various weighting systems are a domain-speific attempt to adapt the raw edit distance metric into something more "perceptual". The needs & error tolerances for one domain, say large vocab open speech, may be different from the command-and-control, "assistant" systems.

Humans differentially weigh different types of ASR errors and are less forgiving with errors (like say insertions) which may distort the meaning of the transcription.

Have a look at this paper: Predicting Human Perceived Accuracy of ASR Systems


Finding amplitudes from FFT data of a 2D image I'm working on experiments with 2D - wrinkled surfaces.
Basically I am getting a bunch of images like this.
This example is hand-made.

Edit: The test image was made by sampling a cosine function something like
$4*\cos{2\pi(3x+4y)}$

Is there a method to find the amplitudes of the first few dominant sinusoids in the image?

I've seen for simple 1-D examples how to get to the amplitudes. It seems to require that the coefficients of the expansion cancel out just right.

Here is a sample image from experiment:

If you can adjust the orientation of your image so that the sinusoid pattern runs vertically to the image, then you can obtain a vertical slice, perhaps at the middle of the image and use that as a 1D waveform for further analysis.

If the waves can come out of your experiment in any orientation you basically have two options:

1) Estimate the angle of rotation ($a$), apply a rotation by $-a$ and then select a slice of the image (just as in the introduction). To estimate the angle of rotation you can look at the output of the Hough Transform which will give you strong "spots" at a given angle for each line (or line-looking) form within your image. That would be the elongated folds of the wrinkles in your image. For more information please see: http://en.wikipedia.org/wiki/Hough_transform

2) Do a straightforward 2D Fourier Transform to recover both orientation and amplitude spectra (at least for the example you provide). The sinusoidal-like form will appear as a linear feature of some orientation in your FFT spectrum which you could pick up with a simple thresholding operation. For more information please see http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf

Hope this helps.


Why notches in sinc filter response Below response shows sinc1,sinc2 and sinc3 filter response.
I understand that notches in the frequency response happens at the output data rate of the sigma-delta ADC.

But I am not able to figure out why there is a notch ? All articles I read mentions that notch happens at integral multiples of output data rate but doesn't give explanation.
Also why notch width increases as order of the sinc filter increases ?

Please let me know if I am missing basic concepts that I should read before finding answer to this question.

Look at the equation for the frequency magnitude response of a sinc filter. It will be of the form:
$$\frac{\sin(\alpha)}{ \alpha}$$
Look carefully at the $\alpha$ term and determine for what frequencies the $\sin(\alpha)$ numerator will be equal to zero. The magnitude response notches occur at those frequencies.


How to perform a gaussian blur using FFT? I am trying to implement a gaussian filter with matlab.Here is my implementation:
So far All I get as output is  a black image.
Any hints?

clear all;
close all;
I = double(imread('Put your path here'))/255;
I=rgb2gray(I);
sigma=1;
[M,N,s]=size(I);
f1=-fix(M/2):ceil(M/2)-1;
f2=-fix(N/2):ceil(N/2)-1;
[fx,fy]=meshgrid(f1,f2);
x=exp(-2*pi*pi*sigma*sigma*(fx/M.^2+fy/N.^2));
If=fft2(I);
If=fftshift(If);
If=If(:,:,1)*x;
If=ifftshift(If)
I=ifft2(If);
imshow(I)

Here is the working code,

clear all; close all; clc;

I = double(imread('Cameraman.tif'))/255; % Divide only if image value in [0-255] range
figure,imshow(I);title('Cameraman original')
%I=rgb2gray(I);  % NOTE: This is only for RGB images, commented otherwise
sigma=3;
[M,N,s]=size(I);
f1=-fix(M/2):ceil(M/2)-1;
f2=-fix(N/2):ceil(N/2)-1;
[fx,fy]=meshgrid(f1,f2);

X=exp(-2*pi*pi*sigma*sigma*((fx/M).^2+(fy/N).^2));
If=fft2(I);
If=fftshift(If);
If=If.*X;
If=ifftshift(If);
I=real(ifft2(If));
figure,imshow(I);title(['filtered image with \sigma =', num2str(sigma)])


What do certain terms mean in EEG processing describing frequency bands? I've recently started a project about detecting driver drowsiness with an EEG and came across some terms I didn't quite understood:

  ...increase in theta activity occur... - ...low voltage theta activity
  during drowsiness induces alpha activity [1]

  ...by theta oscillations [2]

  ...drop of alpha waves

  ...automatically detect so-called alpha spindles [3]

I don't quite get what these terms actually mean and how to recognize them in an EEG signal. Do they talk about amplitude changes? Is there any source describing those EEG signal patterns?

Thanks in advance!

[1] Daniel L. Schacter, EEG theta waves and psychological phenomena: A review and analysis, Biological Psychology, Volume 5, Issue 1, March 1977, Pages 47-82, ISSN 0301-0511, http://dx.doi.org/10.1016/0301-0511(77)90028-X.

[2] Wolfgang Klimesch, EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis, Brain Research Reviews, Volume 29, Issues 2–3, April 1999, Pages 169-195, ISSN 0165-0173, http://dx.doi.org/10.1016/S0165-0173(98)00056-3.

[3] Sonnleitner, Andreas, et al. "Alpha spindles as neurophysiological correlates indicating attentional shift in a simulated driving task." International journal of psychophysiology 83.1 (2012): 110-118.

The short answer: The brain is a huge set of neurons that produce electrical activity as part of their function. We can sense this electrical activity and correlate it with brain states (alert, asleep, awake, meditating, etc). The passing of the brain between states (e.g. alert to asleep) is marked by modifications to the patterns of electrical activity.

The not so short answer: If you are "serious" about processing EEGs, get yourself a copy of Niedermeyer and Da Silva's: Electroencephalography: Basic principles, clinical applications and related fields. It is a monumental book on EEG and will answer all of your questions in excruciating detail. But, you have to give the book the time it needs. Especially sections 2,3,4,5 on the basics of signal generation and acquisition and then 11,12,18,26,30,31 and 56 which talk about modifications caused by certain diseases and during sleep.

Now, the brain is a lump of neurons. Different types of neurons that are there for different types of tasks. Some of them respond to electrical activity (from the nervous system) by propagating or inhibiting electrical activity. We can see this in the lab, by observing how does one single neuron behaves but to date it has been impossible to observe what EACH neuron is doing in-vivo in a functioning brain. This is the "promise" of Neuroimaging.

So, when it comes to electroencephalography, an electrical activity (of the brain) is sensed by an instrument (the EEG aparratus). Actually, in this case, a ridiculously small number of electrodes (a few tenths) is used to sense the combined functioning of a ridiculously large number of neurons (a few billions).

Because of this "imbalance" we see patterns in the electroencephalography measurements that, macroscopically, correspond to brain function. We don't really know what is happening at the neuron level, but we can observe waves of electrical activity spreading throughout the surface of the brain and we hypothesize that this is probably because of the brain's properties of segregation and specialisation (and this, go through the introduction, at least once).

In other words, certain parts of the brain are devoted to certain functions and they "connect" together to exchange information as the brain functions.

Because of this, certain "states" of cognition can be cross referenced with certain electrical patterns.

For example:

If you ask an average someone to close their eyes and relax, their brain will kick into "alpha" wave activity. If it doesn't, then this might suggest that this average someone, is not so average after all.
Theta waves are associated with hippocampal activity. If that activity is found to deviate from what is considered "normal", then this is a sign for things starting to go wrong (please see section 3.2 but also 3.1 for some general remarks.). Theta waves are also associated with another "big mystery" of brain function called "sleep".
Spindles are associated with sleep activity and are correlated with phases of sleep. Lack of spindles, is usually not a good sign for brain function.
As far as "...drop of voltage/waves", "...induces.." are concerned: The waves denote brain states. As you emerge from sleep, the brain makes a transition. So, it is expected that you will stop exhibiting theta waves and sleep spindles and you will come back to having alpha waves and ultimately normal dynamic brain activity. Anything that mentions drop of this wave and increase on that wave, etc, denotes the passing of the brain through different states.

As a side note to this, you probably are aware of various youtube videos claiming to relax, boost creativity, put you to sleep and other claims. I am not sure which and how many of them work but their "basis" is that you can synchronise the brain with an external signal source and induce some sort of brain wave to it. So, this is like saying "well, if the brain produces theta waves when it goes to sleep, why don't we replay some sound source, modulated at theta wave frequencies, to the brain to get it to synchronise and go to sleep?". In fact, certain practices such as meditation are like conscious attempts of a brain's "user" to put it into certain regenerative states. But which of these things work, inferred scientifically, is currently an open question.

Hope this helps.


digital low pass Butterworth filter A digital low pass Butterworth filter that has been designed using Bi-linear transformation has been a pole at $z=0.6$. It is also known that the filter's attenuate (at digital frequency) $\omega = 1.2$ is about $44$ dB. Find the filter order. Give at least one other pole of the digital filter (in Z domain).

I'll try to give you some hints to get you started. First of all, you should know that the bilinear transform is given by

$$s=k\frac{z-1}{z+1}\tag{1}$$

If the analog prototype filter is normalized such that its cut-off frequency is $\Omega_c=1$, then the constant $k$ in (1) is given by

$$k=\frac{1}{\tan\left(\frac{\omega_0}{2}\right)}\tag{2}$$

where $\omega_0$ is the desired cut-off frequency of the discrete-time filter. Some more background on this is given in this answer.

You should also know that the poles of the normalized analog filter lie on a circle with radius $1$ centered at $s=0$ (of course the poles only lie on the left half of the circle). If the filter order is odd as in your example (why?), there must be a pole at $s=-1$. Now you can try to figure out how that pole is transformed to the $z$-plane, i.e.

$$\frac{1}{1+s}{\huge|}_{s=k\frac{z-1}{z+1}}=\ldots\tag{3}$$

From (3), knowing that that pole is transformed to a pole at $z=0.6$, you can determine the constant $k$, and from that, via Eq. (2), you can compute $\omega_0$. Now you know that the filter's attenuation is $3\text{ dB}$ at $\omega_0$, and it is $44\text{ dB}$ at $\omega=1.2$.

The last thing you should know is that with a Butterworth lowpass filter of order $N$, you get approximately $6N\text{ dB}$ attenuation per octave, at least in the frequency range considered here. Now figure out how many octaves there are between $\omega_0$ and $\omega=1.2$, figure out how many dB's difference in attenuation there must be between those two frequencies, and from this compute an estimate of the filter order $N$.


Impulse response of a system from the input-output relation I need to find out the recurrence relation and the impulse response of the following system.

The relation seems to be $$y[n] = \frac{2}{16} x[n] + \frac{3}{16} x[n-1] + \frac{6}{16} x[n-2] + \frac{3}{16} x[n-3] + \frac{2}{16}x[n-4]$$

Since there is no recursion the impulse response $h[n]$ would contain only 5 values (below) being the coefficients?

[2/16, 3/16, 6/16, 3/16, 2/16]

It looks like you answered your own question but to be complete..The impulse response $H(z)=\frac{Y(z)}{X(z)}$ can be found from the difference equation which you have found from the block diagram. So taking the z-transform of both sides of the difference equation we have: $Y(z)=X(z)(\frac{2}{16}+\frac{3}{16}z^{-1}+\frac{6}{16}z^{-2}+\frac{3}{16}z^{-3}+\frac{2}{16}z^{-4})$ and so we have $H(z)=\frac{2}{16}+\frac{3}{16}z^{-1}+\frac{6}{16}z^{-2}+\frac{3}{16}z^{-3}+\frac{2}{16}z^{-4}$.

Now taking the inverse z-transform you get the impulse response to be $h[n]=\frac{2}{16}+\frac{3}{16}\delta[n-1]+\frac{6}{16}\delta[n-2]+\frac{3}{16}\delta[n-3]+\frac{2}{16}\delta[n-4]$, which has the coefficients that you listed.


Fractional/Interpolating Delay Line still sounding glitchy I'm trying to implement a simple digital delay line. I want it so that when the user changes the delay amount, it sounds smooth and not glitchy. Currently I'm implementing a fractional delay line but I'm still getting a very glitchy sound when I change the delay amount.

What else do I need to implement besides a fractional implementation?
I can post code if need be.

EDIT:
Here's some code-

- (void)setDelay:(Float32)dly {
    // get delay amount
    delay = dly*srate;
    if (delay != 0) delay -= 1;

    // Get Delay difference
    Float32 outPointer = w - delay;
    // wrap write
    while (outPointer < 0) outPointer += DELAY_MAX;

    // Integer conversion
    r = (int)outPointer;
    // Wrap read
    if (r == DELAY_MAX) r = 0;
    // Fractional Calculation
    alpha = outPointer - r;
    // Get interpolation values
    omAlpha = 1.f-alpha;
}

// interpolation method
- (Float32)nextOut {
    // If we can get the next output
    if (doNextOut) {
        // First 1/2 of interpolation
        nextOutput = cBuf[r]*omAlpha;
        // Second 1/2 of interpolation
        if (r+1 < DELAY_MAX) {
            nextOutput += (cBuf[r+1]*alpha);
        }
        else {
            nextOutput += (cBuf[0]*alpha);
        }
        // set false
        doNextOut = false;
     }
     // return interpolation
     return nextOutput;
}

If you change the delay parameter suddenly, there will be a jump to another point in time in the data, resulting in a glitch. The easiest is to filter the delay parameter with a one-pole low-pass filter:

$$\text{out}[k] = \text{out}[k-1] + a(\text{in}[k] - \text{out}[k-1]),$$

where $a$ is a constant between 0 and 1, $\text{in}[k]$ is the delay parameter controlled by the user and $\text{out}[k]$ is the delay parameter actually used in the delay line, for integer time index $k$. ($\text{out}[k-1]$ is the previous filtered value.) You can cascade these stages in series to get further smoothing. Two stages should sound good because it will get rid of discontinuities in the pitch shift due to a kind of doppler effect that changing the delay length creates.

The constant $a$ can be tuned manually. A good value might be the highest value that does not produce disturbing glitches. For a better explanation of $a$, useful if you need to adapt to a different sampling frequency:

$$a \approx (1-p)^{\frac{f_s}{t_p}},$$

where $p$ between 0 and 1 represents the percentage that the filtered delay parameter has come towards a new value of the input delay parameter in time $t_p$, and $f_s$ is the sampling frequency. The formula is for a single filtering stage. If you have multiple stages the meaning of $t_p$ changes but you can still use it to go from one sampling frequency to another. The -3 dB cutoff frequency $f_\text{cutoff}$ for a stage will be approximated by:

$$f_\text{cutoff} \approx \frac{-f_s\ln a}{2\pi} $$


Hidden zero in system equation $H(z)$ An FIR linear phase filter has unit sample response $h[n]$ that is real with $h[n]=0$ for $n < 0$ and $n > 7$.  If $h[0]=1$ and the system function has a zero at $z=0.4e^{j\pi/3}$ and a zero at $z=3$, what is $H(z)$?

Solution says, that it has $7$ zeros and here's the system equation:

$$
H(z)=A(1-0.4z^{-1}+0.16z^{-2})(0.16-0.4z^{-1}+z^{-2})\cdot \\(1-3z^{-1})(1-\frac{1}{3}z^{-1} )
$$

My question is:  How does this system function $H(z)$ have $7$ zeros?

When I count them I see:

2 zeros for: $(1-0.4z^{-1}+0.16z^{-2})$

2 zeros for: $(0.16-0.4z^{-1}+z^{-2})$

2 zeros for: $(1-3z^{-1})(1-\frac{1}{3}z^{-1} )$

6 zeros total.

The book says it has 7 zeros. however, i'm not seeing the 7th zero. How to find it?

If the problem description is correct, i.e., if $h[n]$ is zero for $n<0$ and $n>7$, and if we assume that $h[7]\neq 0$ (and why should we assume otherwise?), then the filter has $8$ taps and it is a $7^{th}$ order FIR filter with $7$ zeros. Since it is a linear phase filter and the number of taps is even, it must be either a type II filter (even number of taps, even symmetry), or a type IV filter (even number of taps, odd symmetry). In the first case it must have an additional zero at $z=-1$, in the second case it must have a zero at $z=1$. Refer to this answer for more details on the $4$ types of linear phase FIR filters.

So if we assume even symmetry (type II), the zeros are $z_0=0.4e^{j\pi/3}$, $z_1=z_0^*$, $z_2=1/z_0$, $z_3=1/z_0^*$, $z_4=3$, $z_5=1/z_4$, $z_6=-1$.

If we assume odd symmetry (type IV), the zeros are $z_0=0.4e^{j\pi/3}$, $z_1=z_0^*$, $z_2=1/z_0$, $z_3=1/z_0^*$, $z_4=3$, $z_5=1/z_4$, $z_6=1$.

The corresponding impulse responses are

$$h_{II}=[1,   -5.23,   11.84,  -12.42,  -12.42,   11.84,   -5.23,    1]$$

$$h_{IV}=[1,   -7.23,   24.31,  -48.57,   48.57,  -24.31,    7.23,   -1]$$

EDIT:
The answer above is correct (as far as I can tell) for the given problem description. Now that I've checked the original source (Schaum's Outlines of DSP), which includes the solution, I believe that there are the following possibilities: either there's a typo in the problem description and they actually meant that $h[n]$ is zero for $n<0$ and $n\ge 7$ (note the "greater or equal" sign). In this case we really have a $6^{th}$ order FIR filter and the given solution in Schaum's Outline is correct. The other option would be that they don't know what they're talking about (and I hope and believe that this is not the case).


What version of sinc function to use ? (With additional factor of $\pi$ or not.) I have an exam tomorrow on Signal Processing and one of the questions can go two ways:

Use a Hamming window to design an Nth
th order (N+1 tap) band pass FIR filter with a pass-band
between $f_0$-$f_1$ and a sampling frequency of $f_s$Hz
Following the windowing method of FIR filter design, use a Hamming window to design a
Nth order (N+1-tap) FIR lowpass filter with a cut-off frequency of $f_0$ Hz and a sampling rate of $f_s$Hz.

Now I have a slight problem with how the questions are answered in Mark Scheme for previous paper. For the band pass he uses this definition of $sinc$:
$$sinc(x)=\frac{sin(\pi x)}{\pi x}$$

While for the low pass he uses
$$sinc(x)=\frac{sin(x)}{x}$$

Here are the extracts:

Even in the first case he states that

But then his answers for $h[n]$ follow the $\pi x$ version of the sinc. I have a suspicion that he used matlab to evaluate the numbers and matlab in fact is using the $\frac{sin(\pi x)}{\pi x}$. During the lectures the only way we were taught to define sinc was in the "normal version" without $\pi$.

But then again looking at Wiki article about sinc. It states that the version with $\pi$ is commonly used in signal processing. Could someone tell me what to use and when? Thanks.

Some time ago, on our sister site math.SE, I wrote about this issue regarding the two different definitions of the sinc function.  Here is a slightly modified version of that answer for the benefit of those who prefer to just stay on dsp.SE.

  In the engineering literature, those who define the Fourier transform as
  $$X(\omega) = \int_{-\infty}^{\infty} x(t)\exp(-j\omega t)\ \mathrm dt$$
  tend to use the unnormalized version $\frac{\sin(x)}{x}$, while those who define the Fourier
  transform as
  $$X(f) = \int_{-\infty}^{\infty} x(t)\exp(-j2\pi f t)\ \mathrm dt$$
  tend to use the normalized version $\frac{\sin(\pi x)}{\pi x}$.  My personal preference is for the
  normalized version because the zeroes of the function are the
  nonzero integers, but, like most people, I have learned to live
  with both definitions and figure out which one an author is using
  even if it is not explicitly stated. There are, of course, zealots
  who say that people who use the convention they do not happen
  to prefer are in a state of sin (or did I mean to say a state of sinc?)

  Incidentally, I would like to say that my preference for the
  definition of the sinc (or sine cardinal) function is
  $$\mathrm{sinc}(x) = \begin{cases}\frac{\sin(\pi x)}{\pi x}, & x \neq 0,\\
1, & x = 0,\end{cases}$$
  and not simply $\mathrm{sinc}(x) = \sin(\pi x)/(\pi x)$ the way the OP
  and Wikipedia states it.

I do disagree with Fat32's didactic assertion that if sinc is not defined in a paper or book or on an exam or homework, then one should assume that it is the normalized version that is meant and proceed accordingly. In my experience, it is more commonly the case that if the frequency variable is $\omega$ (radians per second), then the unnormalized version is used while if the frequency variable is $f$ (Hertz), then the normalized version is used, but this is by no means a hard and fast rule that everyone follows.


Opencv calcHist and calcBackProject in plain english i've started experimenting with c++ and opencv because of i want to learn image processing.

Now, my first exercise is to create a skin detector with calcHist and calcBackProject.

But i don't understand few things:

the statistic interpretation of back project, and why it is named "back project"
i have quite good comprehension of what ranges parameter is in calcBackProject. But i'm really stucked with ranges parameter in calcHist function.
For more detection precision i think could be a good thing use backgrojection in more levels: in each channel of r-g-b ans in each channel of h-s-v. But i don't know how i can combine the different results of calcBackProject of separate channels of rgb and hsv.

And i think that my not-well comprehension is caused by lack in theory of what i'm doing with those 2 methods (see the first point). So please explain me in plain english.

See What is Back Projection in the openCV tutorials

An image histogram measures the distribution of colour (and brightness) of the pixels in an image.

If you take an image and identify a region of interest eg. a hand, and calculate the histogram of the pixels in that object.

Then take that histogram and a second image and essentially reverse the process - you pick the pixels in the second image that match the histogram from the first. It's this reverse process that gives it the name back-projection.

You then make the assumption that areas of the image in the second image that have the same colour distribution as an object in the first image are an image of the same (or similar) object.


Autocorrelation of signal with offset If I have a stationary signal $x(t)$ with zero mean and with an auto correlation $r_{xx}(\tau)$, then what is the auto correlation of $y(t)=x(t)+b$ ?

My calculations so far:

$$\begin{align}
r_{yy}(\tau) &= E\{ [x(t+\tau) + b] [x(t) + b] \} \\
&= E\{ x(t+\tau)x(t) + b x(t+\tau) + b x(t) + b^2 \} \\
&=r_{xx}(\tau) + b^2 \\
\end{align}$$

The problem is that the autocorrelation of a constant is a triangular pulse and not just a squared offset. What am I missing here?

Let me just write an answer because what I wrote above about the expectation is not correct and extremely confusing because I didn't show how the $b^2$ drops out.

By definition:

$$\operatorname{cov}(X, Y) \triangleq E(XY) - E(X)E(Y)$$

So, $$\begin{align}
\operatorname{cov}\big((X_{t+\tau} + b),(X{_t} + b)\big)
 &= E\big((X_{t+\tau} + b)(X_{t} + b)\big) - E\big(X_{t+\tau} + b\big)\times E\big(X_{t} + b\big) \\
 &= E\big(X_{t+\tau} \times X_{t}\big) + b^2 - \Big(E\big(X_{t+\tau}\big) \times E\big(X_{t}) + b^2\Big) \\
 &= E\big(X_{t+\tau} \times X_{t}\big)  - \Big(E\big(X_{t+\tau}\big) \times E\big(X_{t}\big)\Big) \\
 &= \operatorname{cov}\big(X_{t+\tau}, X_{t}\big) \\
 &= \gamma(\tau) \\
\end{align}$$

Apologies for confusion.


Rearranging exp terms I have some problems rearranging exp terms to receive a cos or sin. I guess the solution is silly but I‘m not able to spot my mistake.

I know the solution (since i looked it up with wolfram alpha) has to be a $cos(\pi /4n)$ and therefore  $\frac{e^{-j\pi/4n}+ e^{j\pi/4n}}{2}$.
My goal is to rearrange the two terms $e^{j\pi/4n}+ e^{j7/4\pi n}$ (devided by 2) to receive this cos.
One can clearly see that by removing the $\pi/4$ e term on both sites of
$$e^{j\pi/4n}+ e^{j7/4\pi n} =e^{-j\pi/4n}+ e^{j\pi/4n}$$
this remains
$$ e^{j7/4\pi n}= e^{-j\pi/4n}$$ wich is true, when its drawn on the unit circle.

If one doesnt see this trick, my usual approach is to split an e term up to make the remaining two e terms conj. compex of each other. Like this
$$ e^{j\pi n} \cdot (e^{-j3\pi/4n}+ e^{j3\pi/4n})$$

Saddly this results (as far as I can see) in $(-1)^n\cdot 2cos(3/4\pi n)$, which is not what I‘ve expected.
How comes this rearrangement doesn‘t result in the same cos? Is this rearrangement in some cases prohibited?

Hint:

So you have proved that
$$\exp\left({j\frac \pi4 n}\right)+ \exp\left({j\frac{7\pi} 4 n}\right) =\exp\left({-j\frac \pi4 n}\right)+ \exp\left({j\frac \pi4 n}\right)$$

The second derivation does results in the same answer as well. It is not prohibited but one would want to simplify things here. It does not help factorizing with a factor $\exp({j\pi n})$ because $\exp({j\pi n}) = \left(-1\right)^n$ and not $-1$. But one could factor by $\exp({j2\pi n})$ for instance since $\exp({j2\pi n}) = 1$. With this in mind:

$$\exp\left({j\frac\pi4n}\right)+ \exp\left({j\frac{7\pi}{4}n}\right)={\exp\left({j2\pi n}\right)}\cdot\left[\underbrace{\exp\left(-j\frac{7\pi}{4}n\right)}_{\exp\left(j\frac{\pi}{4}n\right)}+ \exp\left(-j\frac{1\pi}{4}n\right)\right]$$
And the two forms are identical.


Is there a term called "Fourier Beamforming"? I understand delay and sum beamforming, where we want to add the signals obtained from different sensors. But adding different time samples from different element of the array depending on the angle we are interested. For broadside or zero degree we would add the samples from the same time instant from each element of the array. On the other hand if we are interested to beamform in a particular angle there would be time delay between the samples of each element.

The question where "Fourier beamforming" comes into this picture? I have seen this term in many places. Any comment would be appreciated.

I'm guessing it's used in the following context.

Assuming you have a physical source with some of an amplitude distribution in a plane: you can calculate the polar pattern of that source simply as the Fourier Transform of the amplitude distribution as a function of space (not time!).

Fourier Beamforming would be the inverse process: you start with the desired polar pattern and do an in inverse Fourier Transform: that gives you the desired in amplitudes as a function of space, i.e. it tells where to put sources or sensors and what weights to apply.


How to find peak value of an analog signal efficiently after sampling in the digital domain? I have a bandlimited analog signal for which I want to find the peak value in real time. The signal is sampled and processed digitally at just enough sampling rate. Since the peak of the analog signal may not be sampled exactly, taking the maximum of the digital samples would give a lower than actual value. So how do I find the peak?
Theoretically, the only way I can think of, is to upsample the signal using an ideal sinc interpolation filter to very high sampling rate, so that we obtain a sample close to the actual peak of the analog signal. Then we can just take the maximum of the digital samples. But this seems roundabout and computationally expensive. Is there is a simpler method to do this?

You are unlikely to be able to evaluate correctly an ideal sinc interpolation filter in real-time (or even with a small time lag), since the future samples of your data are unknown. Plus, when you received your last sample $x[n]$, you generally don't know whether you are on  a peak.  So, it somehow boils down to

how real-time you want to be;
how fast you can estimated you are passed a peak.

There are many local possible interpolations. You can also think about  fast option, and close to real-time, only require to collect the next sample $x[n+1]$, and if knowing $x[n-1]$ you consider that a peak is located between  $x[n-1]$ and $x[n+1]$ (and there values are positive), you can use the standard frequency peak  estimators, which only use the above three samples. For instance:
$$n_{\textrm{peak}} = n+\delta$$
with
$$\delta = \frac{(x_{n+1}−x_{n-1})}{(4x_{n} − 2x_{n-1} − 2x_{n+1})}$$
A primer is provided by famous people here: Fast, Accurate Frequency Estimators, Eric Jacobsen and
Peter Kootsookos, IEEE signal processing magazine, May 2007. The above equation is only one of those this paper provides.


Convolve a minimum phase wavelet with delta peak let's assume I have two time series: one with a delta peak at position $x$ and another time series with a minimum phase wavelet that starts at time $y$ with a maximum peak at position $z$.

Now, when I convolve both time series, is the beginning of my minimum-phase wavelet at position $x$ or is the maximum amplitude of the minimum-phase wavelet at position $x$ ?

I know, that when I convolve a zero-phase wavelet, the maximum amplitude is at the position of the delta function.

I am puzzled..........

If I understand the question properly. Consider the case when the delta is at time zero. In this case the convolution does not do anything to the other signal i.e. it starts at time y and the maximum peak will be at time z.

Now, if we delay the delta to time x, then the output of the convolution will also be delayed by x. So the signal will start at time x+y, and the maximum peak will be at time x+z.

If I haven't understood correctly, you may need to include some figures to more clearly explain the situation.


How would I fit a distribution to this image noise? I have collected some noise data from a dimly lit CMOS image sensor.  The distribution of pixel values is tallied below:-

I'd like to be able to simulate this sensor noise.  How would I fit a statistical distribution to this bell curvish graph? And which distribution? Bear in mind that the distribution is discrete.  From the raw samples, I get mean = 18 and standard deviation = 7 (both approximations).

Since you say your distribution is discrete, if your samples are independent, you could try a multinomial distribution as a model
https://en.m.wikipedia.org/wiki/Multinomial_distribution

You can use either maximum likelihood or a Bayesian estimator.

You could also take the mean and standard deviation and use a Gaussian model with rounding.  A goodness of fit test

https://en.m.wikipedia.org/wiki/Goodness_of_fit

would tell you how good the model is.

If your samples are independent noise generation is straightforward .  A Bayesian model is different because you use another (prior) noise generator as an input to the fixed parameter noise generator.

If your samples are correlated you can use a Gibbs sampler.

You should perform an independence test on your data, of which a Google search will provide a number of candidates

You also might consider boot-strapping

https://en.wikipedia.org/wiki/Bootstrapping_(statistics)

It is also very likely that a physics based noise model exists for your device.


Detecting blobs that connects to any other blob, maybe with OpenCV In the image, I connects to universe, but II and III not.

I need to detect both II and III, also I if possible.

Is it possible with current computer vision libraries?

Flood fill  the background, then hough transform to pick out rough circles II and III


Where is the mother wavelet defined in the Fast Wavelet Transform? Referring to the Fast Wavelet Transform, this transform is implemented as a QMF filter bank. This algorithm consists of high/low pass filtering and subsampling. However, a wavelet transform is typically defined by the mother wavelet - the master basis function that is shifted and scaled.
How exactly does this fast wavelet transform define the mother wavelet? I don't see any mentioning of it.

TL;DR: the wavelet appears at the end of the synthesis filter bank, iterated infinitely.
Theoretically founded, practical and fast DSP tools are often derived from continuous theory: think about how the DFT is derived from the continuous Fourier transform, by discretizing both in time (like the Discrete-time Fourier transform) and frequency  (like Fourier series), which are dual variables.
For a function $\psi(t)$ to be a wavelet, very few conditions (admissibility) are required.  However, discretizing the complementary shift and scale parameters $(a,b)$ of $\psi\left(\frac{t-b}{a}\right)$ is more complicated if one want to remain invertible. And even more complicated to obtain a discrete orthogonal basis.
In other words: apart from a handful of well-known cases (Haar, Shannon, Meyer wavelets), if one have a closed-form formula of a nice admissible continuous wavelet, it is VERY unlikely that it can be discretized as a discrete orthogonal wavelet.
The filter bank structure does the job the other way around. If we assume an orthogonal multiresolution framework, we see that an iterated bank of carefully-designed filters can produce a wavelet analysis. Which wavelet? Theoretically, from a two-scale equation, which generally has no closed-form solution.
You can obtain a good approximation of the wavelet shape with the following procedure. Pick a level $L$. In the deepest detail subband (near the approximation), put a one, and zero elsewhere. Then, do the inverse wavelet transform. The highest the level, the better the approximation.Here is an example with a very q&d code.

dwtmode per
nSample = 1024;
data = zeros(nSample,1);
waveletName  = 'db3';
for iLevel = 2:7
[C,L]=wavedec(data,iLevel,waveletName);
C(L(1)+(L(2)+L(3))/2) = 1;
% C((L(1))/2) = 1;
waveletAtLevel = waverec(C,L,waveletName);
subplot(2,3,iLevel-1)
plot(waveletAtLevel(find(waveletAtLevel)));axis tight;grid on ;
xlabel(['At level ',num2str(iLevel)])
end


Finding the fundamental frequency of a periodic signal Suppose we have the signal $$x(t) = e^{j\omega_1 t} + e^{j\omega_2 t} + e^{j\omega_3 t},$$ where all the frequencies are rationally related (that is, the ratio of any pair of frequencies is a rational number).

How do I prove the fact that the fundamental frequency of the expression above is given by $\omega_0 = \gcd(\omega_1, \omega_2, \omega_3)$?

My attempt is as follows: assume that the signal has a fundamental period $T_0$ which is related to the fundamental frequency by $\omega_0 = 2\pi/T_0$. Now for the signal to be periodic with this period as the fundamental period the following equations must be satisfied $$\omega_1 T_0 = 2 \pi l$$ $$\omega_2 T_0 = 2\pi m$$ $$\omega_3 T_0 = 2\pi n,$$ where $l, m, n \in \mathbb{Z}$ and their values are such that $T_0$ is minimized.

Where do I go from here?
How do I show that the fundamental frequency $\omega_0 = 2\pi/T_0$ is the greatest common divisor of $\omega_1$, $\omega_2$, and $\omega_3$?

1) You are almost there if you have found an $ l,m,n$ and a $T_0$ that fit your equations.  Let $p$ be the largest common factor of $ l,m,n$.  If $p=1$ you are done, otherwise divide $ l,m,n$ and $T_0$ by $p$ and now you are done.

2) From your own equations:

$$ \omega_1 = \omega_0 \cdot l $$
$$ \omega_2 = \omega_0 \cdot m $$
$$ \omega_3 = \omega_0 \cdot n $$

You already know that $l,m,n$ don't have a common factor greater than one.

Hope this helps.

Ced

=================================================

Followup:

Taking it from the top.  You are supposing that you know the frequency ratios:

$$ \frac{\omega_1}{\omega_2} = \frac{a}{b} $$
$$ \frac{\omega_2}{\omega_3} = \frac{c}{d} $$

So $a,b,c,d$ are known.  Therefore:

$$ \frac{l}{m} = \frac{a}{b} $$
$$ \frac{m}{n} = \frac{c}{d} $$

Cross multiply:

$$ lb = ma $$
$$ md = nc $$

Let $ m = bc $ and therefore $ l = ac $ and $ n = bd $

You now have an initial $l,m,n$ set that you can possibly reduce as stated above.  Then the fundamental frequency, $ \omega_0 $, can be solved for using any of the three values.

$$ \omega_0\ = \omega_1 / l = \omega_2 / m = \omega_3 / n $$

Finally, $T_0 = 2\pi / \omega_0 $ gives you the period.

Quite Easily Done

Note, it is possible that none of $l,m,n$ are equal to one.


Difference between recursive and system with memory What's the difference between a recursive system and a system with memory?

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

Is a recursive and with memory, can a system have no memory and recursive?

you can have memory of past outputs.  thems would be $y[n-1], y[n-2]...$  and you can have memory of past inputs which would be $x[n-1], x[n-2]...$  and you know your current input $x[n]$.

if the output you're trying to compute, $y[n]$ will depend on any of the past outputs, it is recursive.  if not, it is not recursive.

if the output, $y[n]$, depends on any past outputs or past inputs, then it is not a "memoryless" system.


Buffering block and filter length in real-time processing I have a signal sampled at 256 Hz, which I want to filter with a 50 taps long FIR filter in real time.
Would it be a problem, if my data block size is only 32 samples?
And should I then concatenate 3 blocks, convolve them with the filter and output only the middle part (to avoid discontinuities)?

And more general questions: were can I read more about all real-time aspects of time-series filtering? In particular about block processing?

Many thanks for help.

Your data block can be as short as 1 sample (although longer is usually more efficient), as long as you save each full convolution response vector to overlap-add/sum into all the following output blocks.  Look up overlap-add/save fast convolution see how this can be done using an FFT/IFFT plus some add/save buffering.


Where should I set my anti-aliasing filter corner frequency for this signal? I have a Fourier transform of an arbitrary digital signal. The transform has For 1024 frequency bins. The magnitude plot of the first 100 is below, with no surprises in the unplotted 924 bins.

Now, suppose that I want to design a single-pole anti-aliasing filter for this signal. How would one approach selecting a corner frequency for that filter?

Suppose if your sampling rate is fs, then the corner frequency of the Anti Aliasing filter would be fs/2


Convolution by IIR filter, a case where circular convolution is allowed? The question is, how to attain the frequency domain representation of any signal (which is finite, non-unit length) convolved with an IIR filter.

The correct answer in my opinion is that there is in fact no answer, which might seem very strange at first. But if we work through the math and use the property of circular convolution we of course will still get an answer.

So let's say our signal is simple discrete step function of the length N-1, starting at t=0, the frequency domain representation is then:

$$F(\omega)=\large e^{j \omega\frac {N-1}{2}}$$

To make the point clearer let's say the IIR filter is of the form. $h(n) = \cos(\omega_0 n)$, which converter to frequency domain is equal to:

$$
H(\omega)=
\begin{cases}
1, & \omega = \omega_0
\\ 0, & \text{otherwise}
\end{cases}
$$

Calculating the frequency representation of the convolution can be done using the formula $$Y(\omega)=F(\omega) H(\omega)$$

Where $Y(\omega)$ is the frequency domain representation of the convolution. However, there is a problem, convolving these functions necessarily results in time domain aliasing, as there is no possible way of having zeroes near the edges. And so, the aliasing occurs, thus the result is:

$$
F(\omega)  G(\omega)=
\begin{cases}
\large e^{j \omega\frac {N-1}{2}}, & \omega = \omega_0
\\ 0, & \text{otherwise}
\end{cases}
$$

In time domain the result is either sinusoidal wave or nothingness, depending on the frequency of the sine wave and width of the pulse we started with.

So continuing this line of thought, sine waves in LTI systems have linear response because of circular convolution (as without the fold back the sine wave would distort). This means that LTI systems properties are a completely theoretical construct, which we of course already knew since infinite sine waves with no end or beginning don't exist in real world. Either way, my line of thought is that you can't convolve an infinite signal using the equation, as no amount of infinities makes a signal fold back onto itself, as using the math will do (admittedly, I may be completely wrong on this). I will be trying to solve the equation for a case where the sine wave is casual next...

A frequency representation of a particular interval can of course be attained (which conviniently cuts the impulse response to a FIR, removing the funny things caused by infinities). The problem also disappears if the filter settles to zero after infinite time, which probably satisfies the definition of IIR, but is of no interest for the purpose of the discussion.

Am I correct on this or does letting the aliasing occur lead to the correct answer for the convolution functions in cases of infinite time functions? If that is the case, what is the rationale behind it?

To answer your first question formulated in the first sentence: if $X(\omega)$ is the Fourier transform of the input signal, and $H(\omega)$ is the frequency response of the IIR filter, the Fourier transform ("frequency domain representation" as you call it) of the output signal is simply $Y(\omega)=X(\omega)H(\omega)$. This is valid in continuous time as well as in discrete time if you use the discrete-time Fourier transform (DTFT) (and not the DFT/FFT).

In the (discrete) time domain, you generally have an infinite convolution sum:

$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]\tag{1}$$

where $x[n]$ is the input sequence, and $h[n]$ is the impulse response of the filter. If either $x[n]$ or $h[n]$ (or both) are of finite length, the sum in $(1)$ becomes a finite sum.

If you manage to clarify the rest of your question, I might be able to add more relevant information to this answer.

If I may guess what confuses you, it may be that you think in terms of the DFT, which can only handle sequences of finite length (and may cause circular convolution artefacts). However, you can implement an IIR filter directly (and exactly, if we ignore quantization effects for the moment) in the time domain. In the general case, the DTFT of the output signal can only be computed approximately. However, for finite length input signals fed into an IIR filter, we can obtain an exact expression for the DTFT of the input signal as well as for the frequency response of the IIR filter. The Fourier transform of the output signal is then simply the multiplication of the two.


frequency spectrum of the conjugate of $x(t)$ Is the frequency spectrum of $x^*(t)$ the same as the one of $x(t)$.

I'm telling this because I found that $$X^\ast(f) = (X(-f))^\ast$$

So to find $X^\ast(f)$, intuitively i would just flip $X(f)$ because of the minus sign and reflip it because of the conjugate. Then flipping two times gives me the original spectrum.

Tell me if I'm doing something wrong.

Thank you very much.

If $X(f)$ is the Fourier transform of $x(t)$, then the Fourier transform of $x^*(t)$ is given by $X^*(-f)$. That's why for real signals $X(f)=X^*(-f)$ holds.


Stability of passive LTI system Is a passive LTI system will always be a stable(bibo and lyapunov) system?
Or in other words
Passivity of a LTI system implies stability?
If above statements are true then there must be some mathematical relation (description) should be there between passivity and stability, so anyone can help me about its mathematical description (if it exists)
And if passivity doesn't implies stability then there must be some counter example to it to prove my assumption wrong.
So whether my assumption is wrong or there exists a mathematical relation between passivity and stability?
I asked this question here
https://electronics.stackexchange.com/questions/513052/stability-of-passive-networks
but I didn't find any useful answer so I posted it here

A counter-example would be an ideal LC-circuit, which is passive but not BIBO stable: if excited at its resonance frequency, the output grows without bounds. Such systems are sometimes called marginally stable, because they have poles on the stability boundary (the imaginary axis).
Note that if a system is not lossless, i.e., if there is a positive resistance such as in an RLC-circuit, the system is strictly passive, and as a consequence it is also BIBO-stable. I.e., strict passivity implies BIBO-stability (all transfer function poles are in the left half-plane).


Why is sum of values of an edge detection filter zero? Usually the filters for detecting edges have the sum of all the values in the filter matrix equal to 0.
For example: [-1 -2 -1, 0 0 0, 1 2 1].

Why is it so ?

An edge detection filter is, by definition, a high-pass filter.  It is looking for quick changes (i.e. high frequencies), not slow trends.  Thus, a good edge detection filter will have a response of 0 at DC.  A FIR filter whose taps sum to 0 has a response of 0 at DC.


Guidance/Brainstorming for a mapping/classification problem I'm quite new to machine learning/data mining and I'm struggling to find the correct path for my problem and would appreciate some guidance or criticism of my proposed solution i.e. is there a better/simpler algorithm for the problem?

The Problem

I have a number of features that describe a particular type (label) of wave (frame of audio) at a predetermined level 'v'. I want to be able to identify the level of an unknown wave and distinguish it from other types of waves that fall under the same higher level category.

Assumptions

A group in a test set should be in increasing order of level v
The type of wave in the group should be the same and known

Proposed Solution

Stage one: Level Selection

For a given type of wave compute the features at each level for N number of samples
For each level calculate the mean/median of each feature the N samples to create a feature vector for each level.
Normalise the feature set by subtracting the empirical mean and dividing by the variance.
Take the Euclidean/Manhattan distance of an incoming vector with the feature set and choose the closest level.
For a group with assigned levels, compare levels with neighbours and report negative differences (should be ascending)  or large differences.

Stage two: type selection

Take the Euclidean/Manhattan distance of an incoming vector with the
feature set for each type at a specific level or maybe across all levels choose closest type.

Extension of Problem

Features evolve over time as well as level

Proposed Solution

Repeat the stages of the above solution for each frame.

Thanks for any help

*Update
I cannot guarantee that the levels v are equivalent across the data I can only guarantee that the order is increasing. i.e. Sample A may have 5 levels v= 1,..,5 and they correspond to {1,..,5} and sample B has 10 samples v = 1,2,..10  and they correspond to {.5,1,1.5,...,10}. How can I capture this without knowing the relationship between levels and identify those which do not follow this pattern. Pleas let me know if this is not clear

You need to be certain that your feature vector is capturing the properties you are looking for. Using the mean and median for dimension reduction (i.e. extracting the features) is not a good choice unless your problem is really simple. Better choices are: PCA, ICA, and SVM's. Capturing the features in the frequency domain can also be a great tool if your signal is sparse in the FFT or DCT domain.


How Can PCA Be Used in Image Analysis I am still a not how PCA can be used in image analysis and where is it is mostly used.  For example how can PCA be used in order to differentiate between different faces? Can you please mention other uses it has?

Imagine you have a set of 10,000 images (32 x 32) of faces.

An intuitive way is to think they have a lot in common.
One step farther would be that if you take one of the faces you could generate it from a Linear Combination of the other 9,999 images.
Sounds reasonable, isn't it?

Let's treat the images as a vector space with dimension $ d = 32 \times 32 = 1024 $.
So all we need is $ 1024 $ independent vectors in order to be able to span any images of size $ 32 x 32 $ in general and faces in specific.

But, if $ 1024 $ images can span all $ 32 x 32 $ images in the world, can't we do better for a subset of this amazing number of different images?
After all we're only after images of faces, human faces.

So we want smaller number of vectors which spans face images.
Maybe not perfectly, but at least almost perfectly based on much smaller number than 1024.
Can we do that?

Yes!
Actually many good ways to do so.
The specific way which uses Principal Component Analysis (PCA) is called Eiganface (At least one of them).
This method is the most efficient way to minimize the $ {L}_{2} $ norm of the reconstruction error when there is a limitation of using $ k $ elements of the linear reconstruction.

If one will analyze the $ {L}_{2} $ norm of the error one would see the best elements are the Eigenvectors of the data set.
Those are the result of the PCA process (Or SVD more generally) and in the context of this method, unsurprisingly, they are called eigen faces.

This general problem is called dictionary learning.
Given a set of data we're learning a dictionary (Set of Atoms) which describe it most efficiently according to some measure of the reconstruction quality.

PCA is one of the most used methods for this kind problems.
The other side of the coin of this problems is the Dimensionality Reduction method where PCA is also heavily used.


Parabolic fiting for sub-sample precision I am performing cross-correlation between two double[1024] arrays. I want to track the correlation's maximum with sub-sample accuracy, so I decided to implement parabolic fitting to my code. I used this class, that seems to work well.

To make this class work I have to feed it with the points I want it to apply to. My question is simple : what are the best points to feed it ? Should I give it the maximum of the cross-correlation and it's direct neighbours ? Should I space them more ? Should I NOT feed it the maximum of the cross-correlation ? Should I limite myself to 3 points, of give it as much as I can ?

Bonus question : if you think this class is not the most suitable and know a better way to implement parabolic fitting to my code, I am open to suggestions.

Thanks.

You can use the formulas presented in the answers to: How to calculate a delay (correlation peak) between two signals with a precision smaller than the sampling period?

To recap, find the largest value, called $\beta$. Take also the values of the samples just to the left of it, $\alpha$, and just to the right of it, $\gamma$. Then calculate the peak position $p$ of Lagrange parabolic interpolation by:

  $$p = \frac{1}{2} \frac{\alpha - \gamma}{\alpha - 2\beta + \gamma}$$

If your signals are white you can get an unbiased estimate $d$ of the peak position from the biased estimate $p$ by:

  $$d = \frac{\sqrt{32p^2 + 1} - 1}{8p}$$

These estimates give the peak position relative to the position of the sample with the highest value.


Oscillations on Hilbert envelope I am trying to extract the decay properties of an oscillating time history. The obvious method for doing this is to use a Hilbert transform. The figure shows an example of a theoretical decay and the envelope calculated using the Hilbert transform.

My problem is that the envelope has oscillations (Gibbs oscillations) due to the fact that there is a discontinuity in the sudden start of the data.

The extracted decay rate, below, which should be a constant in this theoretical case, is unacceptably corrupted by the oscillations.

I have considered filtering the envelope with a low pass filter but as the oscillations are the same order as the original frequency I will spoil the decay data. Alternatively I could create a mirror image of the data and thus remove the starting step. However, this requires finding the starting phase of the time history. Could these methods work or can you suggest another method?

Edit
As suggested by Olli Niemitalo I have had another look at filtering the envelope signal. The filtered envelope is added here

The filtered envelope signal now has no oscillations but has dropped down from the peaks. This is fine for my application were I need the rate of decay. The next figure compares the decay calculated from the filtered envelope to the unfiltered envelope

This is success! Thanks to Olli.

Thanks

A close alternative to a Hilbert transformer is using two infinite impulse response (IIR) all-pass filters that have a phase difference of approximately 90° over the frequency range of interest. For envelope detection, the two all-pass outputs are used in place of the original signal and its Hilbert transform. Here is one such filter pair with quadrature outputs $\text{out}_0$ and $\text{out}_1$:

$$\text{mid}_0[t] = 0.4794008656\,(\text{in}[t-1] + \text{mid}_0[t-2]) - \text{in}[t-3]\\
\text{mid}_1[t] = 0.8762184935\,(\text{mid}_0[t] + \text{mid}_1[t-2]) - \text{mid}_0[t-2]\\
\text{mid}_2[t] = 0.9765975895\,(\text{mid}_1[t] + \text{mid}_2[t-2]) - \text{mid}_1[t-2]\\
\text{out}_0[t] = 0.9974992559\,(\text{mid}_2[t] + \text{out}_0[t-2]) - \text{mid}_2[t-2]$$
$$\text{mid}_3[t] = 0.1617584983\,(\text{in}[t] + \text{mid}_3[t-2]) - \text{in}[t-2]\\
\text{mid}_4[t] = 0.7330289323\,(\text{mid}_3[t] + \text{mid}_4[t-2]) - \text{mid}_3[t-2]\\
\text{mid}_5[t] = 0.9453497003\,(\text{mid}_4[t] + \text{mid}_5[t-2]) - \text{mid}_4[t-2]\\
\text{out}_1[t] = 0.9905991567\,(\text{mid}_5[t] + \text{out}_1[t-2]) - \text{mid}_5[t-2]$$

The coefficients were optimized for minimum largest error over a 20 Hz to 22.03 kHz pass band at a sampling frequency of 44.1 kHz. For a sinusoidal input the largest possible residual error in the envelope is -44 dB or ±0.6%. The following envelope is detected for a unit impulse at time 0:

The vertical axis is magnitude and the horizontal axis is time in samples. Because the all-pass filters are causal, the envelope will not be disturbed until after any transient. Reversing the exponentially decaying input signal gives quite a nice envelope all the way to the transient:

There is a bit of a time lag, which is why the input signal peaks are higher than the envelope. The time lag will be worse for a frequency close to the design band limits of the all-pass filter pair.

If your bandwidth and error requirements are different, you can optimize your own polyphase allpass filters using the HIIR library by Laurent de Soras.

About filtering the envelope: A basic property of linear time invariant (LTI) systems is that if you filter an exponential function then the filter's output is the same exponential function multiplied by a constant. That does not affect the calculated decay rate. The same applies to your truncated exponential envelope if you only use smoothing filter outputs that are not affected by the end effects, by discarding part of the output or by using an anti-causal smoothing filter. The perspective is that you would use smoothing to reduce noise in general, not to reduce end effects. But if you do use it to reduce end effects, I don't think there is an equivalence between envelope filtering and generating the envelope using an all-pass filter pair, because the first does not preserve signal power.


Order Analysis Signal Processing I'm currently working on vibration detection machine used to detect unbalance from a rotating shaft, rotating from 300-2500 Hz. I have both the tacho signal ( 1 PPR ) and the vibration from transducer.

Now I'm really confused. After the transducer there is sweep filter that uses the tacho signal as center frequency so at 1 kHz the output sine wave is 1 kHz.

How do i extract useful data from this signal ?

I've read for the past week on the internet about order analysis and different methods but im so confused about how to implement it, how fast should the ADC sample rate be for the order tracking?

Should I do Fourier transform over the sine wave?

How to transfer the signal to angular domain?

Everything is just confusing. A little help on order analysis, or how to extract information from rotating shaft would be very appreciated!

Most signal processing starts with time-sampling a signal: you have two signals, the vibration signal, $v$, and the tachometer signal, $t$. The first thing you need to be sure of is that these are sampled synchronously in time. Do you know if that's the case?

The reason order analysis is used in the analysis of rotating machinery is that the various frequencies present in the vibration signal often depend on the rotational speed of the main shaft --- which can vary with time.

Order analysis gets around this by resampling the vibration signal synchronously with rotation angle of the main (reference) shaft.

So what you need to do is to take the time-sampled signals you have and produce a single rotation-angle-sampled signal.

To generate the rotation angle signal, you need to take the tachometer signal, $t[n]$, and find the instantaneous phase of it. The way I'd start this is to find the Hilbert transform of $t$ and form the analytic signal:
$$
a[n] = t[n] + j H\left\{ t[n] \right\}
$$
where $H$ is the Hilbert transform and then the instantaneous phase is:
$$
\theta[n] = {\rm atan2}\left( {\tt Im}\{a[n]\} , {\tt Re}\{a[n]\} \right)
$$

From there you should be able to resample both $a$ and $v$ so that you get each evenly sampled in phase.

Once $v$ is evenly sampled in phase, you can then just use an FFT to do order analysis.

  yes they are sampled synchronously in time, does the hilbert transform still need to be applied ?

If they were not sampled synchronously, then it would be almost impossible to perform order analysis at all.

  What is the purpose of the sweep filter ?

The aim of the sweep filter is to turn the tacho signal into as pure a tone as possible.  If it's not a tone, then the Hilbert transform trick will not work correctly for finding the instantaneous phase.

  do i still need to apply FFT? Or does FFT only works on the broadband signal?

Once the vibration signal is resampled, it will be synchronous with the main shaft.  To get the orders (frequencies, in effect), you need to have some sort of frequency plot. The FFT is the easiest.

The effect of this analysis is illustrated in this Matlab application note. I believe the matlab function rpmordermap effectively computes the short-time Fourier transform (but for "orders" rather than "frequencies").


Evaluate expressions without computing DTFT Let $X(\omega)$ be the DTFT of the sequence $x[n]$ given by:
$$
x[n] = \{4, 2, -1, 5, -3, 1, -2, 4, 2\},\quad\text{for}\quad n \in [-6, 2]
$$
I do want to compute

$X(0)$
$X(\pi)$
$\displaystyle\int_{-\pi}^{\pi} X(\omega) d{\omega}$
$\displaystyle\int_{-\pi}^{\pi}|X(\omega)|^2 d{\omega}$
$\displaystyle\int_{-\pi}^{\pi} \bigg\lvert\frac{dX(\omega)}{d{\omega}}\bigg\rvert^2 d{\omega}$

What I have tried

For $X(0)$, I computed the sum and so it was 12.
For $X(\pi)$ I did $X(\pi) = x[9] = \sum_{n = 0}^{8} x[n]e^{-i\omega n} = (-1)^9 \times 12 = -12$, because 0 corresponds to 0, 2 to $\pi/8$, -1 to $\pi/4$, ..., 2 to $\pi$.
For $\int_{-\pi}^{\pi} X(\omega) d{\omega}$, I tried by Parseval Theorem, so it was $2\pi x[n] x[n]*$, but I don't know if it's correct. Last two I don't have any idea of how to do it.

I also have a similar exercise, but it is with DFT instead of DTFT, but it is pretty the same, isn't it?
Thanks in advance.

For the first two questions you just need to use the definition of the DTFT
$$X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{1}$$
and use $\omega=0$ and $\omega=\pi$, respectively.
For 3. you just need to use the definition of the inverse DTFT:
$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega\tag{2}$$
Finally, for 4. and 5. you need Parseval's theorem. In addition, for 5. you need to express the sequence corresponding to the derivative of $X(\omega)$ in terms of $x[n]$. You should be able to either derive this yourself or look it up in a table.


What are the ripple and attenuation in a set of moving average filters? Ladies, Gentlemen,

Consider following set of moving average filters:

(1) - 100 coefficients equal to 1/100
(2) --- 95 coefficients equal to 1/95
(3) --- 90 coefficients equal to 1/90
(4) --- 85 coefficients equal to 1/85
(5) --- 80 coefficients equal to 1/80
(6) --- 75 coefficients equal to 1/75
(7) --- 70 coefficients equal to 1/70
(8) --- 65 coefficients equal to 1/65
(9) --- 60 coefficients equal to 1/60
(10) -  55 coefficients equal to 1/55

Suppose filters are implemented in sequence. What are ripple and attenuation of the whole set of filters?

Responding to Mr Commentator I define ripple 1 dB, cut off band 1 tone (1/6 of octave) and attenuation as little as possible.

Regards

George Theodosiou: As Matt L. said, your final filter is not a standard filter having, say, a visible passband with some measureable passband ripple. For example, the combined frequency magnitude response of just your first three moving average filters implemented one after the other (what we call "cascaded filters") is the following:

The value of 0.1 on the horizontal axis is equivalent to 0.1 times the filter's input data sample rate measured in Hz.

What we can say is: The combined frequency response of multiple filters implemented one after the other is the product of all the filters' complex-valued frequency responses.


Explain the Adaptive Part of Adaptive Algorithms - Kalman Filter and Least Mean Square / Constant Modulus General questions:

Is the Kalman filter (they have used Unscented Kalman Filter) adaptive or not? Is the Unscented Kalman Filter used in the paper an adaptive algorithm?
Adaptive algorithms such as Constant Modulus and Least Squares are adaptive. Why? What is being adapted ? Based on my understanding, the step size is adapted and the weights as well. But I am not sure for the case of Kalman filters -- is the Kalman Gain getting adapted? Adapted to what?

Adaptive Filters are called "Adaptive" when they can adapt to changes in data.
In the filters you mentioned above, which are part of the Linear Filters family the property means their coefficients are changing over time.

Linear Filters are basically weighing and summing the data.
For instance, given no prior information on data you may want to have exact weight for any data given.

Yet in most Signal Processing use cases we'd like to (Or might want to) give higher weight for each data sample according to its properties such as SNR, How good it fit the model, etc...

In the classic Kalman Filter the coefficients are set according to a model matrix - $ P $.
For Kalman filter this matrix stands for the Covariance of the estimation.
It changes according to 2 other matrices which are properties of our model and data - $ Q $ - The model confidence level matrix, $ R $ the data confidence model (For Gaussian data it is basically the SNR).

Since the algorithm support the cases those are changing over time it is called adaptive.

For instance, in the case of tracking a target using RADAR the algorithm can change the matrix $ R $ according to the SNR of the measurement.

The Recursive Least Squares have similar properties (Though it depends only on one factor, something similar to the matrix $ R $ of the Kalman Filter).
One could even see the RLS as a private case of the Kalman Filter (Or the Kalman Filter as a generalization).

There are many extension for the Kalman Filter for many cases.
The UKF (Unscented Kalman Filter) you mentioned is built to handle Non Linear cases of the Model.
But adaptive concept stays the same, the weights are changing according to properties of the Data and the Model.


Determining if a system is Linear Having $h[n]=u_0[n]+0.8u_0[n-1]+1.6u_0[n-2]$ and $x[n]=(u_1[n]-u_1[n-3])$

The goal is to determine if the system $y[n]=h[n]*x[n]$ is linear.

I know that I would have to test it for Homogenity and Additivity in order to determine if the system is linear or not, but the actual proccess testing it is getting me confused.

My professor did it the following way:

$$y_1[n]=h[n]*x_1[n]$$

$$y_2[n]=h[n]*x_2[n]$$

$$h[n]*(\alpha x_1[n]+\beta x_2[n])$$

$$y[n]=h[n]*\alpha x_1[n]+h[n]*\beta x_2[n]$$

Wasn't this insufficient to prove linearity? Wouldn't we have to test it with the actual values of $x[n]$ and $h[n]$?

To show that a system

$$y[n] = T\{x[n]\} $$

is linear, you have to show that it's additive and homogeneous which is codified in

$$T\{ a ~x_1[n] + b ~ x_2[n] \} = a ~T\{x_1[n]\} + b ~T\{x_2[n] \} $$ for all complex $a,b$ and all $x_1,x_2$.

The property set above is satisfied for a system given by :
$$y[n] = T\{x[n]\} = \sum_{k=-\infty}^{\infty} h[k] x[n-k] = h[n] \star x[n] $$ where $\star$ is the convolution operator and $h[n]$ is some sequence to be defined (incidentally it will be the impulse response of the given LTI system).

Of course you didn't have to prove it, as it's already known that if the input-output relation of a system is given by convolution sum; i.e. $y[n] = h[n] \star x[n]$, then it's by definition an LTI (linear and time-invariant) system. Indeed input output relation of any LTI system is given by the convolution sum...


Definition of Ideal Low pass filter (Time Continous) Actually I got confused about definition of frequency response of Ideal low-pass filters because in some books they mention that phase of H(f) should be linear in pass band and it's value is $ -2\pi t_{0}$

However in other books they mention that phase should be zero in pass band.

So I want to know their differences and where does $ -2\pi t_{0}$ comes from?

An ideal low-pass filter has frequency response $H(f)$ equal to the brickwall function: $$H(f)=\mathbb I[|f|<f_c]$$ where $\mathbb I[s]$ is the indicator function, which is equal to 1 if the statement $s$ is true and 0 otherwise, and $f_c$ is the cutoff frequency. Since this frequency response is real, its phase is zero.

This ideal filter is not stable and non-causal. It also has zero delay; it corresponds to $t_0=0$ in the first definition you list above. Now consider the case where we want an ideal LPF, but with some delay $t_0>0$. If the filter input is a sinusoid of frequency $f_0<f_c$, then the filter output is  delayed a time $$-\frac{-2\pi f_0t_0}{2\pi f_0}=t_0.$$

In other words, the first filter in your question is ideal with a delay $t_0$ that you can measure or specify; the second filter is ideal with delay equal to 0. As pointed out in another answer, this helps to make the filter "more causal" -- it allows you to implement a filter that approximates the ideal.

See here for another explanation of phase and delay.


Residual carrier tracking with PLL Residual carrier systems such as PCM/PM/NRZ and PCM/PM/Bi-phase are common in space applications. Compared to suppressed-carrier systems such as BPSK, the literature sources on residual-carrier systems are pretty old, dating back from the 1970s. The demodulators described in these sources are analogue. I'm having some trouble with GNU Radio SDR implementation of these analog structures, in particular, carrier recovery with PLL.

In residual carrier systems, shown above, there is an unmodulated carrier on the Q-channel (pilot) that is used for carrier recovery. The transmitted signal is
\begin{equation}
\mathrm{s\left ( t \right )} = \mathrm{\sqrt{2P}}\sin \left [ \mathrm{\omega_{0}t + \sum_{i=1}^{M}\theta_{i}d_{i}\left ( t \right )} \right ]
\end{equation}
where M is the number of subcarriers and $\theta$ is the modulation index. The modulation index determines the ratio between the power in the data carrier and the residual carrier. Assuming no subcarriers (M = 1), the PCM/PM signal takes the form:
\begin{equation}
\mathrm{s\left ( t \right )_{M=1}} = \mathrm{\sqrt{2P}\cos \theta_{
1} \sin \omega_{0}t + \sqrt{2P}d_{1}\left ( t \right )\sin\theta_{
1}\cos \omega_{0}t}
\end{equation}
In I/Q format, this becomes:
\begin{equation}
\mathrm{s\left ( t \right )_{M=1}} = \mathrm{I\left ( t \right )\cos \omega_{0}t - Q(t)\sin \omega_{0}t}
\end{equation}
where
\begin{equation}
\label{i}
\mathrm{I\left(t\right)} =  \mathrm{\sqrt{2P}s_{1}\left ( t \right )\sin\theta_{1}}
\end{equation}
and
\begin{equation}
\label{q}
\mathrm{Q\left(t\right)} =  \mathrm{-\sqrt{2P}\cos \theta_{1}}
\end{equation}

Im having a problem with implementing the carrier tracking loop, in particular, I dont know what phase error detector to use. The signal in question is shown above, where the residual carrier can be seen as a peak at DC.

For BPSK(suppressed carrier) Costa’s loop, the phase error detector is pretty straightforward and is given by
\begin{equation}
e(kTs) = Q(kTs)*slice(I(kTs));
\end{equation}
Both I and Q channels are used by the PED, as shown in the diagram above. On the other hand, residual-carrier systems use the Q channel only for carrier phase recovery. I have tried different PEDs on my own without any luck.

float
residualPhaseRecovery_impl::phaseDetectorPM(gr_complex out) const
{
    //return imag(out)*slice(real(out)); //Costa's loop, for BPSK
    return imag(out); // Naive Residual carrier PED. Doesnt work.
}

What would be the appropriate PED in this case?

EDIT
Another receiver structure that could be used (?) is shown below (Rice, Digital Communications: A discrete Approach). All I need is an equation for the "Compute Phase Error" block.

I have been facing this problem of phase synchronization with my PCM/PM/SP-L signals as well, and found a temp solution by filtering away the residual carrier.

For carrier offset correction, there are two gnuradio blocks that work, the FLL Band Edge and PLL Carrier Tracking, but I found that the FLL Band Edge blocks works better together with the built-in Costas Loop block when dealing with PM signals with a residual carrier. This is completely from anecdotal experience though, and I have not done any proper characterization of the blocks performance.

One hacky solution of the phase synchronization problem might be to filter away the carrier of the signal before feeding it into the Costas Loop. This filter step should be after carrier offset correction, before Costas Loop. Since the presence of the residual carrier makes the IQ constellation of a Binary PM signal offset from the origin, filtering away the carrier will center the PM signal about the origin, and the Costas Loop PLL can lock its phase to that of a standard BPSK.

Rough Flow Graph:

PCM/PM RX -> FLL Band Edge -> Highpass Filter -> PFB Clock Sync -> Costas Loop -> BPSK Demodulator


struggling to understand why Fourier basis is orthogonal Studying DSP on my own time on Coursera.

Was given a proof to why the Fourier basis is orthogonal, but I can't figure it out. Here is how it is proof goes.

Consider the Fourier basis
$$
\left\lbrace \mathbf w^{(k)} \right\rbrace_{k=0,...,N−1}
$$

defined as:

$$
\mathbf w^{(k)}_n = e^{−j\frac{2\pi}{N} nk}
$$

Let us compute the inner product, that is

\begin{align}
<\mathbf w^{(k)}, \mathbf w^{(h)} > &= \sum_{n=0}^{N-1} \mathbf w^{*(k)}[n] \mathbf w^{(h)}[n]\\
&= \sum_{n=0}^{N-1} e^{j\frac{2\pi}{N} nk} e^{−j\frac{2\pi}{N} nh}\\
& = \sum_{n=0}^{N-1} e^{-j\frac{2\pi}{N} n (h-k)}\\
&= \begin{cases} N & \text{if}\quad k=h\\[2ex] 0 & \text{otherwise} \end{cases}
\end{align}

I do understand all of the derivation except the last step. How does one arrive at the conclusion that when $k\neq h$ the inner product is $0$. Would appreciate any help with that.

Just use the formula for the geometric series (I use $l=h-k\neq mN$):

$$\sum_{n=0}^{N-1}e^{-j\frac{2\pi}{N}nl}=\frac{1-e^{-j\frac{2\pi}{N}Nl}}{1-e^{-j\frac{2\pi}{N}l}}=\frac{1-e^{-j2\pi l}}{1-e^{-j\frac{2\pi}{N}l}}=\frac{1-1}{1-e^{-j\frac{2\pi}{N}l}}=0,\quad l\neq mN$$


What is the Asymptotic Equipartition Property (AEP)? I am currently studying about Polar Codes in 5G standard and while reading my paper I found something called AEP which is required for channel coding. I surfed the web but didn't found a satisfying answer. Can someone explain what it is, clearly?

When you are sampling a stochastic process $n$ times, the larger you make $n$, the higher the probability that the series of samples is contained in the so called strongly typical set of outcomes of length $n$, where all members have roughly the same probability to be realized. If you let $n\rightarrow \infty$ it holds that
$$-\frac{1}{n}\log p(X_1,X_2,...,X_n)\rightarrow H(X)$$
with $H(X)$ being the entropy rate of the process.

This should be intuitively clear: the longer the sequence, the more possible outcomes, each with possiblity near zero but positive. The number of typical outcomes with $p\approx\epsilon$ is growing faster than the number of non typical outcomes as $n$ is growing. As $n\rightarrow \infty$, the "leftover probability" for non typical outcomes gets smaller and smaller, thus making the above equation valid.


Method to remove small unwanted region and fill holes for specific image I have a binary image that maintains some unwanted region (small white dot) and hole regions (in figure 1). My idea is that to fill the object after I will remove the unwanted region but I don't now exactly wich method can help me

For this image in particular? I believe it should be sufficient to threshold (to get a binary image with only the visible particles present), then perform area opening (remove all objects with an area smaller than a certain value). Then, depending on how the remaining object looks, perform morphological operations such as thickening (simply enlarge the object) and closing (to fill in holes).


determine two signals with a scale factor Suppose I have 2 signals from function $f_1(x)$ and $f_2(x)$, respectively, and assume the sampling rate is above Nyquist frequency, so we can restore the underlying functions $f_1(x)$ and $f_2(x)$. But my question is, from the 2 signals, how to tell if the underlying functions $f_2(x)$ is scaled from $f_1(x)$, i.e. to determine if $f_2(x)=f_1(ax)$, where $a$ is a non-zero real number.

Thanks!

A necessary (but not sufficient) conditions for $f_2$ to be a temporally scaled version of $f_1$ is that a spectral representation with a logarithmic frequency scale (such as the constant-Q transform) of $f_1$ is a translation of a log-frequency spectral representation of $f_2$.

Practically, given two signals, you can perform the test and evaluate $a$ by computing the CQT of $f_1$ and $f_2$, cross-correlating them and looking at the location of the peak. The strength of the peak might give you an idea of the spectral similarity of the two signals irrespectively of their temporal scale; and the position of the peak will give you the temporal scaling factor.

This type of representation robust to temporal scaling is useful in music signal modeling, where the different notes produced by a music instrument are - in a very rough approximation - temporally scaled versions of themselves.


First order Wiener–Hopf filter design Consider a random process with auto-correlation function:
  $$r_{\rm dd} [k] = \beta^{\lvert k \rvert}\quad\text{where}\quad 0 < \beta < 1.
$$
  Suppose also that the observation is:
  $$
x[n] = d[n] + v[n]
$$
  where $v[n]$ is uncorrelated white noise with variance $\sigma^2$.

  Design a first order Wiener–Hopf filter to reduce the noise in $x[n]$ of the form

  $$W(z) = w(0) + w(1)z^{-1}$$

I refer to the notation from the Wikipedia article.

Your received signal is $x[n]$, and its Autocorrelation is given by
$R_x[n]=R_d[n]+R_v[n]$ when noise and signal are uncorrelated. Hence,

$$
R_x[n]=\beta^{|k|}+\sigma^2\delta[n]
$$

The cross-correlation between the received signal and the signal of interest is

$$R_{xd}[n]=E[d[n'](d[n'+n]+v[n'+n])]=R_d[n]$$

under the assumption again that signal and noise are uncorrelated.

Now, the Wiener-Hopf equation gets you

$$
\begin{pmatrix}R_x[0] & R_x[1]\\ R_x[1] &R_x[0]\end{pmatrix}\begin{pmatrix}w[0]\\w[1]\end{pmatrix}=\begin{pmatrix}R_{xd}[0]\\R_{xd}[1]\end{pmatrix}
$$

Filling in the variables, we get

$$
\begin{pmatrix}1+\sigma^2 & \beta\\ \beta &1+\sigma^2\end{pmatrix}\begin{pmatrix}w[0]\\w[1]\end{pmatrix}=\begin{pmatrix}1\\\beta\end{pmatrix}
$$

Now, you just need to solve for $w[0],w[1]$ to get your filter coefficients.


When a stochastic process would be a beneficial model in terms of noise Let's say we have an image/signal with some noise in it. When would it be beneficial to model the signal as an outcome of a stochastic process?
More specifically:

How significant would noise have to be in for example an image for it to be beneficial to estimate the power spectrum density instead of taking the fourier transform directly?

And would it matter what kind of noise that was present in the image/signal, if it is additive or multiplicative? I don't think I'm seeing how a stochastic process could be a beneficial model for a signal yet.

I hope these weren't stupid questions. I'd be really grateful for any answers to them!

When would it be beneficial to model the signal as an outcome of a stochastic process?

When the process that generates the image/signal has a strong element of chance.

This is not related to noise necessarily. Sometimes, a deterministic model for the behaviour of a quantity simply does not exist.

For example, modern supermarkets have barcode readers that go "beep" when customers scan an item. Suppose now that we want to create a model of that "beep" sequence from which we can derive the rate at which customers are bying products from the supermarket.

What is a sequence of "beep"s? It is Beep-Random Pause-Beep-Random Pause-Beep and so on.

There is no deterministic equation for the inter-beep pauses. Deterministic models of the form $y = \alpha t + \beta x \dots$ return the same output $y$ for the same parameters $\alpha, \beta, \dots$.

But in the case of the supermarket, we have no idea about how customers arrive at the barcode scanner to scan their baskets. They just arrive...at random.

So, in this case, the inter-beep random interval is governed by a Poisson distribution which only requires a rate parameter. This rate denotes average rate of events per unit of time. In our case, the rate at which customers scan their products.

As this model is stochastic, it generates a different series of interbeep intervals everytime that you "run it" because there are infinite ways by which you can generate a sequence of events at some average rate $\lambda$...at random.

The exact same mathematics govern radioactive decay, the output of a suitable detector (that detects byproducts of the decay) and after a few more parameters are taken into consideration, the progressive formation of a gamma camera image.

Now, a stochastic process does not have to produce discrete events at some rate $\lambda$ only. It might represent the random fluctuation of some other quantity. What is the birth rate of a society at a given year?.....I don't know....But let's say that it is a random variable with characteristics governed by some probability distribution. Random variables represent quantities that either are inherently uncertain (e.g. radioactive decay) or, our understanding during the construction of the model is uncertain.

With this in mind, you might find exploring the rest of the examples provided in Stanley Pawlukiewicz's post really entertaining too. Drunken sailor walks are everywhere around us and they are the starting point of a lot of wonderful generative art and graphics effects today.

  How significant would noise have to be in for example an image for it to be beneficial to estimate the power spectrum density instead of taking the fourier transform directly?

The choice between obtaining the Power Spectral Density (PSD) or Discrete Fourier Transform (DFT) is not exactly based on the amount of noise present in the image but rather, what are you trying to model.

  And would it matter what kind of noise that was present in the image/signal, if it is additive or multiplicative?

Yes it would. It is the difference between an additive noise component, like background crowd murmur in a shopping centre when you are trying to hear a person talking to you and a fault at the cable of the microphone recording the person. In the first case, you record two sources simultaneously and both are present. In the second case, when one is present the other is not. If your microphone cable is faulty, it kind of makes and does not make contact so your rBZZZord is fuZZZBZ f artifacBZZZSPAK!ZZZ ike thiZZZZ. You either have the person's voice or not. Well, when do I have the person's voice? I don't know...what is the rate at which the microphone connection fails at? Another manifestation of multiplicative noise is wow and flutter in record players. That is, the variation that is inserted in the recording because of the variable rate of rotation of the disc platter. There, the spectrum of the playback is the modulation (or product) of the song's spectrum with the variation of the rate of turn of the disc platter. What could that be? I don't know, but let's say that it's a random variable that is $0$ (that is, bang on the right speed), plus or minus $\sigma$ turns per minute.

In fact, if we were talking about disc players that use an elastic middle wheel to couple the motor to the disc platter, this "model" would directly represent the variation in the radius of the elastic wheel that progressively gets worn out and inserts this "modulation" of the player's speed.

I am mentioning this last part more to show that there is a direct link between reality and how you choose to represent it in a model.

So, the real question here is "What are you trying to model"? And one of the best places to look for experience in creating models of the real world is Automatic Control.

Hope this helps.


What is Bottom-Hat transformation? I am reading the Digital Image Processing book by Gonzalez and i can't seem to understand Bottom-hat transformation. Can anyone explain it to me in easy terms?

The explanation in Wikipedia for Top-hat transform - Properties seems pretty simple, for the bottom-hat, another name for the black top-hat.

The black top-hat returns an image, containing the "objects" or
"elements" that:

are "smaller" than the structuring element,
are darker than their surroundings.


Is there a way to evaluate block size when doing image processing or video processing I am currently doing foreground segmentation in image sequences. A DCT 8 by 8 feature is extracted for each pixel.

Now i want to explore dynamic block size for different scenes. Is there any rules or criterion to decide optimal block size?

Thank you!

I assume the block size is related to image quality no matter you are doing the segmentation or compression. Generally speaking, smaller blocks are better as far as correlation (the pixels in the block are more correlated), but if they're too small, you won't be able to eliminate many coefficients from each block without causing severe distortion, especially when there are high variance and more details in the block. Some compression algorithms use different block sizes within the same image so that high-detail regions are divided into smaller blocks, while smooth regions are divided into larger blocks (Chong et al's Compression of medical images through adaptive block-size DCT coding, and The Influence of DCT Block Size on Coding Efficiency: An Adaptive Variable Block Size Transform Coding System with Lagrangian Cost Function Decision Criterion).

Besides, Karunasekera et al's A distortion measure for blocking artifacts in images based on human visual sensitivity provides you a reference to estimate the distortion induced by blocking artifacts. Pramod el al's Unsupervised Segmentation of Medical Images using DCT
Coeﬃcients applied DCT coefficients to the image segmentation.


What is facebook's image compression algorithm? For a project I'm working on (master thesis), I need to figure out what is facebook's image compression algorithm.
The goal: to be able to mimic facebook's image compression locally (using python).
The images I use have the following exif: Image exif (Q = 90)
After facebook compression, I have the following exif : Compressed image exif
Note: If the images I uploaded with the link seem unavailable, please refresh the page! :)
After searching for a while, I could not find a recent and clear answer to my question.
I tried to run the compressed image through ImageMagick's identify. It tells me that it is a JPEG compression (DCT progressive + Huffman) with a quality of 85. But when I try to do my own compression using python, the result is not the same.
Here is my code:
from PIL import Image
img = Image.open('original.jpg')
img.save('output.jpg', optimize=True, quality=85, progressive=True)

Do you have any ideas or resources that could help me? I'm running out of ideas.
Thank you!

If you have access to the bitstream, you should be able to figure out which codec they are using.
Assuming that it is good old jpeg, possibly with a tweaked encoder implementation, I would parse the stream using JPEGSnoop to learn as much as possible. It can even try to match «known»
encoders:
https://www.impulseadventure.com/photo/jpeg-snoop.html
If you can upload uncompressed bmp images to facebook and observe the first-generation compressed jpeg, you might learn even more as you have a reference.
Reverse engineering an encoder black box like this may be tedious unless they went with a plain vanilla implementation.
JPEG is actually a flexible standard, and «quality=85» is a simplified way of expressing how 8x8 DCT coefficients of luma (and chroma) are quantified. Sniffing the stream will allow you to read out those 2x64 quantization scalings, and also (if any) chroma subsampling, and (unlikely) non-bt601 color primary matrixing. An encoder is allowed to set any (small) dct coeff to zero if it likes (in the interest of smaller files), besides the regular quantization. Also, the encoder may use non-standard huffman tables but that should affect the bitstream not pixels.


What's the maximum allowable amplitude modulation frequency for a given carrier frequency? I'm currently studying wavelets and had an interesting thought experiment: If you were to calculate the wavelet transform of a signal using a wavelet of a fixed frequency, you would get time varying wavelet coefficient. Say this fixed frequency wavelet is 1000Hz. You essentially yield a signal that shows how a 1000Hz frequency component varies with time.

If you were to run a Fourier transform on this new signal, what is the maximum frequency content you could find in it?

My intuition says that there is some upper limit to this (i.e. you couldn't find 10kHz frequency content), but I'm curious if there is a way to prove this limit.

After thinking about it some more, it seems like what I'm asking is: what the maximum modulation frequency that can be encoder into a carrier frequency?

In the frequency domain, an AM modulated sinusoid looks like a carrier plus two side-bands in complex conjugate symmetry above and below the carrier.  If you permit the lower sideband to cross into the negative frequency spectrum, then there is no limit to maximum modulating frequency.  Of course that may not look like your typical AM modulated signal in the time domain.


Image transformation where a curve becomes a straight line Suppose that I have the green curve (image attached for visualization), values stored in a vector and I want to transform the image such that the curve becomes a straight horizontal line. So if the curve becomes a straight line all normals of the curve has to be parallel to the vertical axis.

How can I achieve such a transformation?

I would recommend a grid transform. Each point of the grid has a $(dx,dy)$ offset, and a point in the middle of a square is interpolated according to the control points. Specifically in your case, you can keep most of the points in place, except those which are near the curve you want to straighten.


extending windowing function to higher dimensions I was looking at this particular thread that discusses extending windowing functions for singling function to higher dimension. There is one bit in the accepted answer that I could not get my head around.

So, they create a multi-dimensional windowing function but then they scale the weights in this line.

# scale the window intensities to maintain image intensity
np.power(window, (1.0/data.ndim), output=window)

What is the purpose of doing that?

Without being familiar with the library they use, this pretty much looks to me like a numerical normalizing condition.

Why you need this

This is needed in order to not change the magnitude of the numbers in your data set after filtering. Consider for example a large (multidimensional) region in your image data that exhibits the same value everywhere, e.g. 529. If you filter/smooth this dataset with a filter that has a FWHM that is considerably smaller than the size of the region, you would expect that the value in this region does not change. If you wouldn't include this normalizing condition, this behavior is not guaranteed and you could end up with totally different values... e.g. 1E6 or 0.02 instead of 529. This would seriously complicate further computation, if you had a threshold or so in later processing steps. Or if you converted your data later on... The values would always depend on the filter settings.

A somewhat analytical example

Consider a 1D Gaussian curve. It has a normalization factor $\displaystyle\frac{1}{\sqrt{2\pi}\sigma}$ that is needed to satisfy $\displaystyle\int G(x)\,\textrm{d}x = 1$. For higher dimensions this becomes a little bit more complex: $\displaystyle\frac{1}{\sqrt{(2\pi)^k |\Sigma|}}$, with $|\Sigma|$ being the determinant of the covariance matrix.

A somewhat numerical example

If you do not know the analytical formula, you could also just construct your multidimensional Gaussian without the prefactor, numerically integrate it, and use the resulting factor as the normalization constant.


Inputs for system identification What effects does have a different input such as impulse, step, sine sweep etc. on performing a system identification?

I mean what is up/down side, and which one is best, and why chose one rather than the other.

For now I see the sine sweep as the most reliable for determining a system since the other ones doesn't seem to hit nearly all modes of a system.

But what do you say? I am not an expert at all and I am a bit confused here.

Both Impulse and Step Response can help identify a system, however it might be hard to actually perform the identification. You'd have to fourier transform the response. Basically, this is only possible using digital signal processing. Using a sine sweep, the result is easily read, but the sweep will take a relatively long time.


Unit delay operator as parameter in function I'm busy trying to implement a filter from this paper. I'm a bit confused with one particular step (30e in the paper).

Basically I need to find the zero frequency gain (K) of the following whitening/comb filter:

The relevant steps in the algorithm are:

and

where $q^{-1}$ is the unit delay operator, $\omega_0$ is the fundamental freq of the comb and $\rho$ is a constant.

I'm not sure how to calculate K. Since the function isn't acting on anything, what should be delayed? The signal (at time t)? Or does passing a constant in place of the delay mean it should be replaced?

Any help is greatly appreciated!

EDIT

Some clarification on $q^{-1}$:

  The right-shift operator $q^{-1}$ acts on a discrete-time signal f(k, T) as $q^{-1}$ f(k, T) = (k − 1, T). The operation $q^{-1}$ therefore shifts the discrete-time signal to the right by one sampling interval.

With $A(q^{-1},t)$ as given in the formula you get

$$A(1,t)=\prod_{k=1}^n(1-2\cos(k\hat\omega_0(t))+1)\tag{1}$$

and

$$A(\rho,t)=\prod_{k=1}^n(1-2\cos(k\hat\omega_0(t))\rho+\rho^2)\tag{2}$$

So $K(t)$ is simply the quotient of (1) and (2).


Computations in Grayscale morphology dilation I red this link, but unfortunately i still don't understand the process...

By the way, Dilation process is performed by laying the structuring element H on the image I and sliding it across the image in a manner similar to convolution.

The grayscale morphological dilation formula is written as follows :

$[I ⊕ H](u,v) = \displaystyle\max_{(i,j)∈H}\{I(u-i,v-j)+H(i,j)\}$

If we assume a greyscale image I of size 2x3, such that I=[1 2 3; 3 7 2]; and a structuring element H of size 1x4, such that H= [0 0 7 7].

What are u,v and i,j? Please can someone explain to me in details how the formula is computed at each element of I in order to obtain the new image dilated_I ? Because i still don't understand well the concept..

What was meant by a structuring element that contains only the origin ? i didn't understand well the concept behind grayscale dilation.
Any help will be very appreciated.

Structuring element is mostly supposed to be a binary array (though you can try 0 0 7 7 by yourself). H=[0 0 1 1]
u, v are the pixel coordinate at I. In your example, the only pixel value that will change is I(1,2) (the array label starts from 0 instead of 1).

I =  1     2     3
     3     7     2

Consider I(1,2), from the formula, I(1,2) I(1,1), and I(1,0) will be considered.

I(1,2) = max( I(1-0,2-0) + H(0,0), I(1-0,2-1) + H(0,1), ...
              I(1-0,2-2) + H(0,2), )
       = 7

That's the pixel value at (1,2) after dilation.

I will show you I(1,1):

I(1,1)=max(I(1-0,1-0)+H(0,0), I(1,0,1-1)+H(0,1))

That is still the original I(1,1).

EDIT

Mathematically, SE can has any values in its elements, but in Matlab, imdilate only accepts binary SE as a mask to determine the dilation region. You can still use [0 0 7 7] in your own calculation of dilation.

As your suggestion, I am showing you an example with SE of 3X3 size. I will calculate I(2,2).

I =  1     2     3
     3     7     2
     3     7     2
H =  1     0     0
     0     1     1
     1     0     1

I(2,2) = max( I(2-0,2-0) + H(0,0), I(2-0,2-1) + H(0,1), I(2-0,2-2) + H(0,2),...
              I(2-1,2-0) + H(1,0), I(2-1,2-1) + H(1,1), I(2-1,2-2) + H(1,2),...
              I(2-2,2-0) + H(2,0), I(2-2,2-1) + H(2,1), I(2-2,2-2) + H(2,2))
       = 7