Fourier Transforms, Convolution, Cross-correlation: what is their physical unit exactly? Let us assume we are talking about real, deterministic, electrical signals $x(t)$ and $y(t)$ (magnitude in Volts).
There are different kind of Fourier Transforms. I made a table to summarize:

NB: By the "<-" symbol, I mean variable substitution.
I tried to figure the physical unit of the output. The $\mathrm{V}\cdot\mathrm{s} = \mathrm{V/Hz}$ for the FT is fine, but I am not satisfied with what I get with the alternative transforms... is that right?
Also, analyzing the units of convolution product of $x(t)$ and $y(t)$:
$$
(x \ast y)(t) = \int\limits_{-\infty}^{+\infty}x(t-u)y(u) \ \mathrm{d}u
$$
... or the cross-correlation of $x(t)$ and $y(t)$:
$$
(x \star y)(d) = \int\limits_{-\infty}^{+\infty}\overline{x(t-d)}y(t) \ \mathrm{d}t
$$
...it would yield $\mathrm{V}^2 \cdot \mathrm{s}$ unit, which is... I don't know... somewhat closer to an energy?
But, as I understand it, the convolution product is actually a signal, typically the output of a filter bank from the original signal.
------------------------------------- EDIT --------------------------------------
In electronics, for the convolution product, the $y$ function should rather be seen as a pattern. Typically  $y$ is an impulse response (often noted $h$), its unit is $s^{-1}$ and then the unit of the convolution product is $V$, which is legit for a signal.
-------------------------------------- END --------------------------------------
On the other hand, the cross-correlation is rather meant to be some kind of inner-products series, so it makes more sense to see it as some kind of energy (the auto-correlation of a signal at $d=0$ multiplied by some coefficient is actually its energy).
------------------------------------- EDIT --------------------------------------
In the cross-correlation product, both $x(t)$ and $y(t)$ are indeed signals.
The above definition is for finite-energy signals, but it changes for finite-power signals:
$$
(x \star y)(d) = \lim_{T\to\infty}\frac{1}{T}\int\limits_{-\frac{T}{2}}^{+\frac{T}{2}}\overline{x(t-d)}y(t) \ \mathrm{d}t
$$
-------------------------------------- END --------------------------------------
So physically speaking, maybe there are some unitary normalization coefficients missing in those formulas?
------------------------------------- EDIT --------------------------------------
For finite-energy deterministic signal $x(t)$, we get:
$$
R_0.E_x = \int\limits_{-\infty}^{+\infty}|x(t)|^2 \ \mathrm{d}t <{+\infty}
$$
For finite-power, we get:
$$
R_0.P_x = \lim_{T\to\infty}\frac{1}{T}\int\limits_{-\frac{T}{2}}^{+\frac{T}{2}}|x(t)|^2 \ \mathrm{d}t <{+\infty}
$$
We may divide those by $R_0 = 1 \Omega$ to get respectively energy or power quantities (hypothetical, relative to 1 Ohm).
Same goes for cross-correlation, depending on whether the signals are finite-energy and/or finite-power.
-------------------------------------- END --------------------------------------

I agree that the units for Fourier transform are V/Hz. Most of my experience with the convolution integral dealt with a convolution of the impulse response of a linear network with the voltage input to the network. The impulse response for, say, a simple RC network will have the product RC in the denominator. Thus, the unit of impulse response is per second. So, the units of a convolution would be volts-seconds * per second = volts. For correlation, either auto or cross-, in the case of power signals (as opposed to energy signals), you should divide the integral by the period, T. If it is not periodic, still divide by T, but take the limit as T goes to infinity. Now the period, T, has units of seconds and you are integrating over time, which has units of seconds. The time in the numerator and denominator cancel, and you are left with volts squared, which gives power when divided by resistance (assuming that the volts are rms).


The benefits of a fractionally spaced equalizer I often read that one of the benefits of FSE vs a symbol rate equalizer, besides avoiding any aliasing in the roll off region, is that it is  less sensitive to the timing phase. What does that really mean? Does it mean that if at the output of the matched filter instead of taking one sample (assuming it is not exactly at the center of the eye diagram) , two samples are chosen, then the performance in terms of SNR will improve? How? by interpolation? I don't really understand the relation between the equalizer and what is called the timing phase ( not sure what that expression means either)

The more samples you have per symbol the less timing error you have since delay in samples is converted to delay in time by dividing by the sampling rate. I beleive the far more dominant reason for a fractionally spaced equalizer is the aliasing you mention as can be seen in the frequency response of the equalizer. For example, 2 samples per symbol can capture all the phase and amplitude distortion represented over the channel bandwidth uniquely as depicted in the graphic below. (This is assuming the signal was sufficiently filtered prior to being decimated to the fractional rate used).


Getting rid of Canny outliers I have been using the Canny edge detection function in OpenCV to detect the edge of an elliptic annulus (light object, dark background). In theory there should be two edges. However, it sometimes does not all edges on the outer ring.

So if I radially plot my detected points it looks like this (x: angle, y: radius from center of mass)

There is a wavy line which is the outer border of my object and a few outliers which are the inner border (you can tell because they look like they are systematically offset). I don't know the exact breadth of the annulus.

Q: Is there any way to programmatically get rid of these "obvious" outliers so that the exact outer contour is interpolated?

Please see my answer to this question: Detecting and isolating part of an image

You can use circular shortest path to do this quite elegantly. You are using OpenCV and and not MATLAB so some modification of my code is required.


Comparing similarity of two simultaneously recorded EEGs in freq domain I have two EEG signals that are simultaneously recorded during the night (sleep). Imagine one device recorded from channel F4 and the other from Fp2. Now I would like to check the similarity of these signals, so, the first step is to sync them which has been done successfully using cross-correlation :)
I can also report some similarity metrics in time-domain, e.g. Pearson correlation per window (e.g. each 30 s) to quantify the similarity of signals.
The problem arises when I want to compare (or let's say quantify) them in the freq domain. When I plot e.g. spectrogram or PSD we can see their similarity visually; however, I don't know how to quantify this similarity in the frequency domain. By vision the spectrograms (or PSDs) are quite similar but I guess it is not a powerful reason without reporting any similarity metric as a number. Is there any metric by which I can report similarity of PSDs or spectrograms window-by-window and based on each frequency bin (Delta, Theta, etc)?

You can express correlation in the frequency domain but the closest thing to what you are describing would be Coherence.
It returns a number between 0 and 1 per each frequency bin of the spectrum. The closest this number is to 1 the more "similar" the signals would be and you can average this number over the whole spectrum or over a frequency range.
Hope this helps.


Amplitude cutoff following FFT suppose I have a noisy sine wave (as R code):

# Sampling frequency
fs <- 1000
# Sampling period
ts <- 1/fs
# Length of signal
l <- 1500
# Time vector
t <- seq(0, l - 1)*ts
# Signal
s <- 0.4*sin(2*pi*50*t) + sin(2*pi*120*t)
# Noise
x <- s + 2*rnorm(length(t))
# Plot signal
qplot(t, x) + geom_line()

and want to determine the amplitudes of the signal with:

# FFT
y <- fft(x)
# Two-sided spectrum
p2 <- abs(y/l)
# Convert to one-sided spectrum
p1 <- p2[1:(l/2 + 1)]
p1[2:(length(p1) - 1)] <- 2*p1[2:(length(p1) - 1)]
# Define frequency
f <- fs*seq(0, l/2)/l
# Plot spectrum
qplot(f, p1) + geom_line()

How can I decide whether a peak in the spectrum is indeed the amplitude of the signal and not noise? Is there a formula to define a threshold?

Thanks!

You may have to read on SS (spectrum sensing) techniques. There is a famous method called "Energy detection".

Else

If the noise is WGN, then its PSD is constant. So at a particular frequency of the spectrum, if the PSD is around that of noise, then it's noise, else it's signal.


How to compute dBFS? I asked myself how to compute dBFS (dB full scale) from a value of sample between 1 and -1?

and in general?

It is really quite simple for values between 1 and -1:

valueDBFS = 20*log10(abs(value))

If you have values between another range, for example 16bit, then it's:

valueDBFS = 20*log10(abs(value)/32768)

(because 16bit signed has values between -32768 and +32767, e.g. 2^15=32768)

And so the general formula is:

valueDBFS = 20*log10(abs(value)/maxValue)

The code examples are given in MATLAB code.


replace input signal with δ will have the impulse response of the function? Is my impulse response right?
By definition,the impulse response is the output when the input is a impulse signal,so

$y[n]=\sum\limits ^{n}_{k=-\infty}\frac{1}{2^{n-k}}\ x[k]$,the impulse response of $y[n]=\sum\limits ^{n}_{k=-\infty}\frac{1}{2^{n-k}}\ \delta[k]$ ,Is my thinking right?

Of course, if you have a formula describing the system's response to an arbitrary input $x[n]$, you will obtain the response to an impulse if you choose $x[n]=\delta[n]$.

However, I suppose that the idea of the exercise is to obtain a closed-form expression for the impulse response. You should recognize that the given input-output relation is just a convolution sum, and from that you can just write down the impulse response in a very simple form, without using an infinite sum.

You can also find that simple expression from your second formula by realizing that $\delta[k]$ is only non-zero for exactly one value of $k$. Depending on the value of $n$, that value may or may not be reached within the given summation limits.


Which filter for an audio equalizer I’m new to DSP and I have to do an audio equalizer in c++.
I did a lot of research about it  and tried some stuff in the last month, but I’m a little overwhelmed with all those informations and it’s not working yet.

I decided to start over and for now, all I got is a function which fills the audio buffer with a sin wave (1kHz). And I’m trying to apply a filter with those values to it.

<Filter>
    <Bands num-bands="28">
        <band0    freq="31.25"      width="7.23625"     gain="0" />
        <band1    freq="39.3725"    width="9.1171"      gain="0" />
        <band2    freq="49.6063"    width="11.4868"     gain="0" />
        <band3    freq="62.5"       width="14.4725"     gain="0" />
        <band4    freq="78.7451"    width="18.2342"     gain="0" />
        <band5    freq="99.2126"    width="22.9737"     gain="0" />
        <band6    freq="125"        width="28.945"      gain="0" />
        <band7    freq="157.49"     width="36.4684"     gain="0" />
        <band8    freq="198.425"    width="45.9473"     gain="0" />
        <band9    freq="250"        width="57.89"       gain="0" />
        <band10   freq="314.98"     width="72.9368"     gain="0" />
        <band11   freq="396.85"     width="91.8946"     gain="0" />
        <band12   freq="500"        width="115.78"      gain="0" />
        <band13   freq="629.961"    width="145.874"     gain="0" />
        <band14   freq="793.701"    width="183.789"     gain="0" />
        <band15   freq="1000"       width="231.56"      gain="6" />
        <band16   freq="1259.92"    width="291.747"     gain="0" />
        <band17   freq="1587.4"     width="367.579"     gain="0" />
        <band18   freq="2000"       width="463.12"      gain="0" />
        <band19   freq="2519.84"    width="583.495"     gain="0" />
        <band20   freq="3174.8"     width="735.157"     gain="0" />
        <band21   freq="4000"       width="926.24"      gain="0" />
        <band22   freq="5039.68"    width="1166.99"     gain="0" />
        <band23   freq="6349.6"     width="1470.31"     gain="0" />
        <band24   freq="8000"       width="1852.48"     gain="0" />
        <band25   freq="10079.4"    width="2333.98"     gain="0" />
        <band26   freq="12699.2"    width="2940.63"     gain="0" />
        <band27   freq="16000"      width="3704.96"     gain="0" />
    </Bands>
</Filter>

So here are my questions :

In which domain should I apply the filter ?
I started in the frequency domain, but this blog (here and here) shows me that it’s not the best and only solution.
Which kind of filter should I use ?
I’m lost with all those high pass, low pass, peaking, shelving filters; and I don’t know which one to choose.
I only found “tutorial” to create a filter which modify only one frequency / band. Is it possible to create a filter that is a sum of multiples filters ? Or applying multiples filters in series to the signal ?

Feel free to ask if I forgot something that might help to understand my problem.
Thanks for your time.

PS : Sorry if this post is not very specific, but I want to clear that out before posting specifics problems or code.

Simply put, you need a bank of passband filters.  You feed your signal through each of the filters, and sum up the outputs from the filters.

Designing the filters is where the fun comes in.

First off, assuming this is just audio (music or the like) then there's no need of special filters.  You can use the simplest and fastest and not worry too much about phase and ripple.
You will want to use IIR filters, as they require the least amount of calculation - since you have 28 bands and presumably are working in stereo, you will want to reduce the load on the CPU else the computer won't be able to do much else.  IIR uses less CPU time than FIR filters.

There are online tools that you can use to create the needed coefficients for the filters.  Matlab or Gnu Octave can also be used.
You will need to look up algorithms for implementing IIR Filters, and how to make use of the calculated coefficients.

Specifying the bands is the next fun step.
When a filter passes audio, it doesn't simply remove everything out side of its passband.  It just isn't possible to do so.  Frequencies outside of the passband of a filter sort of slope off down towards zero amplitude.
At some frequency outside of its passband, the audio amplitude will drop by one half.  This is the -3dB point, and is used in specifying filters.
What you want to do is to arrange your passbands such that the -3dB points of adjacant filters are at the same frequency.
Say you have two filters.  One should pass frequencies from 500 to 1000Hz and the next from 1000 to 2000Hz.  You would want to make the high cutoff of the lower filter to be 1000Hz, and the low cutoff of the high filter also be 1000Hz.  Design the filters such that 1000Hz is down by half in each filter (usually, just specifying it as the cutoff is enough.)  When you now sum the outputs of the two, they add up to just the audio from 500 to 200Hz, with everything in that range at (more or less) the same amplitude as before.

Now, with your audio split into banks, you can multiply the output of each filter with some constant to dampen or emphasize that band.
This image demonstrates what I mean about the way the bands have to be layed out.


Voice Audio Detection algorithm I have to detect speech intervals in large pre-recorded files in my project. I think, there won't be so much noise in the background (audio will be recordered in the room or even in the studio), but it still can be.

So, what are good VAD algorithms that I can use (algorithms, not libraries)?

Depending on the actual recordings, the algorithm complexity could range from dead easy to really complex...

I'll take the studio recording case first, so I can assume :
- (Almost) no noise coming from outside (cars, trucks, bus...)
- Nobody slamming the door in the middle of the recording
- Voice samples are recorded at optimal level independently of who is speaking (some ppl speak very low, others loud and clear)

If that's the case, then you could implement the logic of what sound engineers call a noise gate : below a certain audio power threshold, it's silence, otherwise it's not.
You'd have to calculate mean power over a very small window and check it against the threshold. Probably take safety margins around "speeched" areas so you don't cut the start/end of the words.

Determining the correct threshold is a manual process involving trials and errors.

Noise gates sometimes use 2 different thresholds : one to determine speech has started and another to determine speech has stopped. The second one being lower than the first one so it doesn't stop in the middle of a word.

Determining threshold(s) for different files recorded at different levels (because of different speakers) might become tricky. One way to help is to first normalize the volumes : i.e. check for the max value and apply gain in order to bring this value to the max authorized one in the range you're working with (might be -1.0 to 1.0 or 0.0 to 1.0 or signed/unsigned integer values).

Another thing to define is the minimum length of "speeched" areas and the minimum length of silent areas.

If there's a bit more background noise (and I assume the noise volume doesn't change too much over time), you could use the property of white noise when calculating autocorrelation factors : it's mean value goes close to zero

Here is an example of "brute force" autocorrelation

int n = vocalSample.length;
for (int j = 0; j < n; j++) {
    for (int i = 0; i < n; i++) {
        correllation[j] += vocalSample[i] * vocalSample[(n + i - j) % n];
    }
}

then calculate mean value of the array.

Do this for very small arrays (e.g. 1ms at 22KHz -> 22 values), otherwise it will take ages.

I actually used this technique with some success on Ted talks to remove silence in order to perform speaker identification.

If you have hours and hours of recording, you might consider multithreading your calculations :-)
In my calculations (remove noise, normalize, extract vocal print), I crunched 10Gb in about 8 minutes using 8 worker threads.

State-of-the-art vad is something beyond what could be explained here, here is a link to an extremely interesting IEEE article on the subject : http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6403507&punumber%3D97


how to take a second order transfer function without any zeros and represent it using two first order transfer functions does anybody know how to take a second order transfer function with no zeros (no s terms in the numerator):
$$ g(s) = K\frac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2} $$
and represent it using two first order transfer functions.
the layout for my second order transfer function is
     (number)/((number)s^2 + (number)s + (number))

this second order transfer function was previously in the form
(number)/((number) s + (number))^2 before i expanded the brackets

For $|\zeta| \le 1$, let $\zeta= \cos\theta$, so $\theta=\mathrm{arccos}\,\zeta$
$$\begin{align*}g(s) &= K\frac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2}\\
\\
&= K\frac{\omega_n^2}{s^2 + 2\omega_n s \cos\theta+ \omega_n^2(\cos^2\theta +\sin^2\theta)}\\
\\
&= K\frac{\omega_n^2}{(s + \omega_n \cos\theta)^2+ \omega_n^2\sin^2\theta}\\
\\
&= \sqrt{K}\frac{\omega_n}{s + \omega_n \cos\theta + j\omega_n\sin\theta} \cdot \sqrt{K}\frac{\omega_n}{s + \omega_n \cos\theta  -j\omega_n\sin\theta} \\
\\
&= \sqrt{K}\frac{\omega_n}{s + \omega_n e^{j\theta}} \cdot \sqrt{K}\frac{\omega_n}{s + \omega_n e^{-j\theta}} \\
\\
\end{align*}$$
For $|\zeta| \ge 1$, since $\cosh\theta \ge 1$, let $\zeta= \pm\cosh\theta$, so $\theta=\mathrm{acosh}\,\left(\pm\zeta\right)$
$$\begin{align*}g(s) &= K\frac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2}\\
\\
&= K\frac{\omega_n^2}{s^2 \pm 2\omega_n s \cosh\theta+ \omega_n^2(\cosh^2\theta -\sinh^2\theta)}\\
\\
&= K\frac{\omega_n^2}{(s \pm \omega_n \cosh\theta)^2- \omega_n^2\sinh^2\theta}\\
\\
&= \sqrt{K}\frac{\omega_n}{s \pm \omega_n \cosh\theta + \omega_n\sinh\theta} \cdot \sqrt{K}\frac{\omega_n}{s \pm \omega_n \cosh\theta  -\omega_n\sinh\theta} \\
\\
&= \sqrt{K}\frac{\omega_n}{s \pm \omega_n e^{\theta}} \cdot \sqrt{K}\frac{\omega_n}{s \pm \omega_n e^{-\theta}} \\
\\
\end{align*}$$


Express a wave as a sum of cosine I started off with this:

$$ x = 5 \cos(200 \times 2 \pi t ) + 10 \cos( 400 \times 2 \pi t)$$

by performing a 1024 point DFT on it. Then I performed IDFT on the complex values obtained through DFT. Now, I want to express the original wave in terms of it's harmonics, so when I combine the DC and first harmonic or upto n harmonics, I get the original wave back. But how do I obtain these harmonics?

This is what I actually want to do. I have some real valued data(144 points). I want to express this data as a sum of cosines. For this purpose, I perform a 144 point DFT operation on this data, and bring them to frequency domain.
Does the Euler's identity, give me these cosine?

$$X[k] * e^{\tfrac{j2\pi kn}{144}} + X[144−k] * e^{\tfrac{−j2\pi kn}{144}}$$

Please help me obtain the Amplitude and phases of these harmonics.

The Fourier transform converts signals into coefficients of cosines and sines. This is why the coefficients are complex. Simply adding up those two coefficients won't give you sufficient information about the phase (someone can correct me on this statement, but I believe it's true). In other words, for signal of length N you're trying to find coefficients $a_k$ and $b_k$ such that

$$ x[n] = \sum_{k=0}^{N/2} a_k \cos(\omega k)+\sum_{k=0}^{N/2} b_k \sin(\omega k)$$

I would suggest using Discrete Cosine Transform (DCT) instead. It literally breaks you signal up into a sum of cosines, it gives you back all real coefficients and it takes care of the phase problem completely.

If you must use Fourier Transform, you should realize that

$$e^{i \omega t} = \cos (\omega t) + i \sin(\omega t)$$

This tells us that we can recreate the even-symmetric component (cosine) from real coefficients of the Fourier Transform, and odd-symmetric component (sine) from imaginary coefficients of the Fourier Transform.

The formula then befomes:

For a given $k$,

$$a_k \cos\left( \frac{2 \pi kn}{N} \right) = \frac{1}{N}\left( \Re (X[k]) e^{ \tfrac{i 2 \pi kn}{N}} + \Re (X[N-k]) e^{ \tfrac{i 2 \pi (N-k)n}{N} } \right), $$

$$b_k \sin\left( \frac{2 \pi kn}{N} \right) = \frac{1}{Ni}\left( \Im (X[k]) e^{ \tfrac{i 2 \pi kn}{N}} + \Im (X[N-k]) e^{ \tfrac{i 2 \pi (N-k)n}{N} } \right), $$

where $\Re(\cdot)$ and $\Im(\cdot)$ represent taking real and imaginary part of a complex number respectively. Note that this is done for each $n \in \{0 \ldots N-1\}$

Note that because of circular symmetry of DFT, $X[N] = X[0]$, so you'll end up adding up the same thing twice according to the formula above. If you do that, then the actual $a_0$ is half of what you get according to the formula above. If you find this division by 2 confusing, this section of Wikipedia article on Fourier Series may provide a better explanation (look for the formula where $a_0$ is divided by 2).


How does the earphone cable in a mobile phone pick up the FM reception? As I understand, the RF frequency-modulated signal centered in the range of $88-108\textrm{ MHz}$, is picked by the earphone antenna and is directly sampled by the sound card adc at $48\textrm{ kHz}$ and then digitally processed. Will this process not cause aliasing? Could anyone please correct my understanding?

Your question is based on a misunderstanding. The RF signal is not directly sampled by the sound card.

Instead, the earphone cable just doubles as an antenna for a dedicated FM receiver; the soundcard hasn't got anything to do with FM reception.
In fact, it doesn't even "get" radio frequencies – they're filtered out before and fed to the FM receive chain.


Services for phoneme recognition? I wonder if anyone can point me to services or dev libraries that focus on or suit phoneme recognition? Target audience is children at learning to read stage.

We have trialed IBM Watson (not good) and Google (only a little better).

Hope this is a suitable question for this forum. Thanks in advance.

SoapBoxLabs have very good technology for children, including pronunciation scoring


Adding AWGN noise with a correct noise power to the signal I would like to add AWGN noise to my BPSK and QPSK signal. So assume the BPSK signal is x1 and QPSK signal is x2. The EbN0 is $10\textrm{ dB}$. Therefore I will first convert the EbN0 from $\textrm{dB}$ to linear scale. Assume Eb=1 and Es=1.

I would like to know can anyone explain to me, should the noise power be N0 or N0/2?

I am really confused since according to this link
N0 is used to add noise to one dimensional noise and N0/2 is used to add noise to two dimensional signal. Therefore if Eb=1 and EbN0=10 $\textrm{dB}$

Eb=1;
N=1e6;

EbN0=10;
EbN0_lin=10^(EBN0/10);
N0=Eb/EBN0_lin;
noise=sqrt(N0)*randn(1,length(x1))

or should it be

noise=sqrt(N0/2)*randn(1,length(x1))

The same thing for QPSK, is it

Es=1;
EsN0=10;
EsN0_lin=10^(EsN0/10);
N0=Es/EsN0_lin;
noise=sqrt(N0)*randn(1,length(x2))

or

noise=sqrt(N0/2)*randn(1,length(x2))

I was thinking that according to convention the N0 is sigma^2/2 however the link above is saying different things.

Can anyone explain the correct way of calibrating the noise?
And explain if the above reference is incorrect or not?

The kind of AWGN channel you're simulating is sometimes called a "discrete time" channel. That means that you're going to simulate the bit error rate at the output of a decision device, the input to which is produced by a correlator. This is more easily explained by looking into the transmission and reception of a single bit. I'll assume BPSK.

Let's say that you transmit the physical signal $a_0p(t)$, where $a_0$ is either $A$ or $-A$, and $p(t)$ is a pulse of unit energy. To make things easier we'll assume $p(t)$ is symmetrical, that is, $p(t)=p(-t)$ (this is true for all practical systems I know). Note that, in this scenario, the bit energy is $E_b=A^2$.

Assuming that there is no noise, the received signal is $r(t)=a_0p(t)$. The receiver correlates $r(t)$ with the transmitted pulse $p(t)$. The output of the correlator is $$\int_{-\infty}^\infty a_0 p(t)p(t)\,dt=a_0,$$ where I've used the fact that the pulse has unit energy.

Now assume that the input to the correlator is pure Gaussian white noise $n(t)$ with PSD $N_0/2$ and zero mean. Then, the correlator's output is $$n=\int_{-\infty}^\infty n(t)p(t).$$ By the Wiener-Khinchine theorem, $n$ is a zero-mean Gaussian random variable with variance $\sigma_n^2=N_0/2$.

In the case of having signal plus noise at the receiver's input, the output will be $$r=a_0+n$$ because the correlator is linear. Using the value of $r$, a decision can be made on whether $A$ or $-A$ was transmitted.

We define the SNR at the correlator's output as the bit energy divided by the noise variance: $$\text{SNR}=\frac{a_0^2}{\sigma_n^2}=\frac{2E_b}{N_0}=\frac{2A^2}{N_0}.$$

Notes:

When you plot the BER, you'll likely use $E_b/N_0$ in the horizontal axis. The reason is that in quadrature modulations, the total noise power is $N_0$, not $N_0/2$.
Practical systems use a matched filter instead of a correlator. The results are exactly the same.
In a practical system there will be delay, you transmit multiple pulses instead of only one, etc. Those details are not essential to understanding how the SNR is measured in your simulation.


Computationally Efficient Ways of Determining If Pixels Are Clumped I'm working on an application where I get a video frame that's already segmented into good and bad pixels (as in -- the hardware either coughs up a pixel value or a flag that says "bad"). I want to come up with a figure of merit for the frame that's based both on the number of bad pixels and the amount that they are clumped, with one big clump being worse than a bunch of little ones.

Here's two pictures to illustrate: the left picture has 17 dots 10 pixels in diameter -- it's better, on the theory that the dots are artifacts that can be "seen around".  The right picture has about the same number of bad pixels, but they're all in one dot 41 pixels in diameter -- it's worse, on the theory that having so many bad pixels in a clump may occlude something interesting.

Things I can think to do:

Walk the perimeter of each clump (really inefficient)
Low-pass filter (originally I was thinking to use an FFT and look for high frequency content (not terribly efficient, and not terribly accurate, either, I think))
draw a box around the bad pixels by finding the largest vertical and horizontal extents, and compare the area of that box to the bad pixel count (this could have false negatives if there's a large blob and a few little dots)

Are there known good ways to do this?

I've not tried this, but I'm wondering if something like a simple low pass filtering, followed by a column sum will do.

No time to explain code right now, but this seems to do something like what you want. Will return in a few hours and explain what it's doing.

R Code Below

library(imager)

N <- 128

noise <- array(runif(N*N*1*1),c(N,N,1,1)) #5x5 pixels, 1 frames, 1 colours. All noise
small_clumps <- as.cimg(noise)
blurry <- isoblur(small_clumps,5)

layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(colSums(blurry-mean(blurry)))
plot(small_clumps)
plot(blurry)

large_clumps <- small_clumps
large_clumps[65:75, 65:75] <- 1
large_clumps[15:25, 35:45] <- 1
blurry_large <- isoblur(large_clumps,5)
# par(mfrow=c(3,1))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(colSums(blurry_large-mean(blurry_large)))
plot(large_clumps)
plot(blurry_large)


Definitions for amplitude spectrum and phase spectrum I want to help someone with a homework in signal processing, and the courses he provided are some slides, without clear definitions in them. I am required in solving a problem to represent the amplitude spectrum and the phase spectrum of a given signal. I searched the internet and didn't found the definitions of these two. I want to create a Matlab script for these, but if I don't know the definitions I cannot do that. Please give me the definitions for these.

  Amplitude spectrum
  Phase spectrum

If there is a Matlab command, or some other program which plots these spectrums directly, I would be very interested. Thank you.

Beni I think you want something like this:

samplingFrequency = 1000; % / Hz
signalFrequency = 50;     % / Hz
signalLength = 500;
signal = cos(2*pi*signalFrequency./ samplingFrequency .*[0:1:signalLength - 1]);
[hh ff] = freqz(signal, 1, 2048, samplingFrequency);
figure(10);
plot(ff,10*log10(abs(hh)));
ylabel('Amplitude in log scale');
xlabel('Frequency / Hz');
figure(20);
plot(ff, angle(hh));
ylabel('Phase');
xlabel('Frequency / Hz');

You do have MATLAB correct?... if you do try it and let me know how else I can help.

Well, a couple things. First off, remember, since this is a digital signal, it means of course that it was sampled at some frequency $samplingFrequency= 128 or 127$ in your case. The physical frequency of your signal in this case seems to be $signalFrequency = 21$.

The ratio of $\frac{signalFrequency}{samplingFrequency}$ is equal to the $f_0$ in your case. This is called the 'digital frequency'.

For your problem it sounds like he also wants you to fix the fft size, (change the 2048 in my code to 128) and the rest of the numbers as I described and run the code. Does he mention how long he wants your signal to be?


Statistical properties of the Kalman estimates under Gaussian noise For a linear state-space model with independent Gaussian state and output noises and perfect guess for initial state, do Kalman estimates have the following properties:
$$
E(\hat{x}_{k|k} - x_k) = 0
$$
$$
P_{k|k} = Var(\hat{x}_{k|k} - x_k),\text{ or }Var(\hat{x}_{k|k}),\text{ or }Var(x_k) ?
$$
where

$x_k$ is the state at time $k$, which is random
$\hat{x}_{k|k}$ and $P_{k|k}$ are Kalman esitmates, i.e. outputs of Kalman filter.

Are there references mentioning these?

Thanks!

The following two statements are equivalent of saying:

$$
E(\hat{x}_{k|k} - x_k) = 0
$$

(1) That the estimator is unbiased; and

$$
P_{k|k} = Var(\hat{x}_{k|k} - x_k)
$$

(2) That the estimator is consistent.

Both of these conditions are necessary in order for the filter to be optimal - i.e. the best possible estimate of $\mathbf{x}_{k|k}$ with respect to some criteria.

If (1) is not true, then the mean-square error (MSE) would be the bias plus the variance (in the scalar case). Clear, this is larger than the variance only and hence suboptimal.

If (2) is not true (i.e. the filter-calculated covariance is different to the true covariance) then the filter will also be suboptimal.  Since the Kalman Gain is based on the calculated state covariance, an error in the covariance will lead to an error in the gain.  Error in the gain means a suboptimal weighting of the measurements.

(As it happens, both conditions are true for a properly modelled filter. Errors in modelling, such as the dynamic model or noise covariances will also render the filter suboptimal).

Source: Bar-Shalom, especially Section 5.4 on page 232-233.


What's the best criterion for determining a frequency peak? I'm writing an algorithm that analyzes results of discrete Fourier transform (DFT). The algorithm should detect amplitude peaks in predefined ranges of frequencies. This is an example of such a peak:

Vertical red lines are range boundaries (left and right) and range center (middle).

What are the best criteria for determining the absence of this peak in a specified range?

For now I'm using standard deviation and comparing it to some calibration value. But it can be easily confused with some high amplitude noise like this:

First convolve the FFT sequence in the given range by a fairly wide gaussian window. In case there is only one peak, this will just widen it. In case there are many peaks (as in the second example), this will merge them into a uniform and less spiky mass.

Then you can use one of the following peakedness metrics:

Kurtosis.
Ratio of geometric to arithmetic mean


Number of discrete samples required for the longest wavelength I am trying to understand the effect of the critical Nyquist frequency when applying the Goertzel algorithm for estimating the power spectrum of a discrete signal.  (Goertzel doesn't really matter, FFT applies as well).  Assume that I generate $N=512$ discrete values of the sine wave $x[n]$, where $n=1,2,\ldots,N$ using the function

$$x[n]=\sin(30 \cdot 2\pi n/N) + \sin(60 \cdot 2\pi n/N)  + \sin(120 \cdot 2\pi n/N) \\  + \sin(260 \cdot 2\pi n/N) + \sin(330 \cdot 2\pi n/N)$$

After applying the Goertzel algorithm, there will be 512 power spectrum values at $k=1,2,\ldots,512$. Looking at the results, it appears that in the upper half of the scale ($257 \leq k \leq 512$), the real (vs. imag) power spikes for frequencies 1/260 and 1/333 are negative.  In addition, in the lower half of the range ($1 \leq k \leq 256$), the power value for $k=60$ is zero, and there are false positive spikes at spurious frequency values.  However, if I double the length of the generated signal to 1024, and only apply the Goertzel algorithm to $k=1,2,\ldots,512$, all of the power spectrum values at $k=30, 60, 120, 260$, and $330$ are positive and there are no false positive values below 512 and the value at k=60 is non-zero and positive.  My understanding of the Nyquist requirement is that the power value at $k=1$ is unreliable.

So what I believe I have observed is that if you want to use the real (vs imag) power spectrum values at frequencies between 1/2 and 1/512, you need to provide a signal with twice (1024) the number of samples.  Another way of saying this is that if you have a discrete signal with 1024 samples, you can only determine the power spectrum at frequencies greater than 1/512, i.e., 1/511,1/510,...,1/2.

Is there a law or equation that states that the sample size $N$ needs to be at least twice as long as the longest wavelength to be assessed?  Or, is this really the meaning of Nyquist, i.e. $2N$.

If it is a real signal, as you are doing, then the negative spectrum will be the complex conjugate of the positive spectrum. Since there are N samples in frequency given N samples in time, and the frequency axis extends from 0 to N-1 where 0 represents DC and N represents the sampling rate (with the upper half of this from N/2 to N-1 equally representing the negative frequency axis) then yes in this case you need 2N samples to represent the spectrum without aliasing; and yes this is explained by the Nyquist-Shannon sampling theorem: https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

However if you used an equivalent complex signal, then N samples would be sufficient.

$$x[n]=e^{j(30 \cdot 2\pi n/N)} + e^{j(60 \cdot 2\pi n/N)}  + e^{j(120 \cdot 2\pi n/N)} \\  + e^{j(260 \cdot 2\pi n/N)} + e^{j(330 \cdot 2\pi n/N)}$$


What does it mean by "Adding Gaussian Noise to a 2D Image"? As far as I understand,

Suppose, we have a 2D image of pixel resolution 200 x 300. That means, the image has 60000 pixels in it.

Now, we would generate $n$ random values and add those values to $n$ number of random pixels.

Now, my question is, what would make those $n$ values to look like Gaussian Random Values?

What would be the logic that I should use?

You can find Gaussian noise generators in many image processing softwares.At each pixel, you add a realization from such random noise generators. You can look for more details at How to Generate White Gaussian Noise.

On StackExchange, additional sources of information:

Why Gaussian noise is usually used?
Why is Gaussian noise called so?
What's gaussian noise look like?
The model of noise in Image denoising problem


PSD using MATLAB's $\tt pwelch$ algorithm window length and measure units When using pwelch algorithm in MATLAB for a signal with $N$ samples,

Is there a good rule of thumb for determining the length $L$ of the sections (second variable for the function)?
And assuming my input is a voltage signal , what are the measure units I get for $10\log_{10}(P_{xx})$ where $P_{xx}$ is the output signal?

Regarding your questions:

  Is there a good rule of thumb for determining the length $L$ of the
  sections (second variable for the function)?

$L$ (or window argument) specifies a window function to apply, prior to FFT computation of the signal. If $L$ is a vector, that vector is considered as window coefficients and is applied on (multiplied to) signal, for example:

Pxx=pwelch(x,[0 0 1 2 3 3 2 1 0 ]); % Means compute spectrum of x.*[0 0 1 2 3 3 2 1 0]

if it is an integer, then a Ham. or Han. window with length of $L$ is generated and then this Hamming window is applied to the signal:

Pxx=pwelch(x,10); % Means compute spectrum of x.*(10 point hamming window)

pwelch is based on FFT, and since the spectrum (DFT) of segmented signals shows frequency leakage, I think it would depend on how much resolution you you want and how much frequency leakage, you can tolerate, but generally The more longer the window for frequency static signals is the more leakage is reduced.
See this great tutorial for more detail:

http://www.fhnw.ch/technik/ime/publikationen/2012/how-to-use-the-fft-and-matlab2019s-pwelch-function-for-signal-and-noise-simulations-and-measurements

A side note for example codes above is that if the length of your signal x is greater than specified length for the window, then your signal is segmented to vectors with length of $L$

  assuming my input is a voltage signal , what are the measure units I get in the output?

Since your signal is pure voltage (unknown impedance and hence no current (A) is associated), it is not actual "power" or "energy", however you might use $\mathrm{V}^2/\mathrm{H}z$ for PSD.
This link might shed more light on it:
https://en.wikipedia.org/wiki/Spectral_density#Electrical_engineering
In case of logarithmic PSD, you've mentioned it is generally regarded as $\mathrm{dB}/\mathrm{Hz}$.


Routh's stability criterion: zeros of the auxiliary polynomial and of its derivative When using RH criterion and using the auxiliary equation special case, which of the following is true?

The auxiliary equation $A(s)=0$ gives some(or all) of the symmetrical poles.
The differentiated auxiliary equation $\displaystyle \frac{\mathrm{d}A(s)}{\mathrm{d}s}=0$ gives some(or all) of the symmetrical poles.

I have seen both of them mentioned either separately or together in different books. I am really confused on which one is actually correct.

After solving a few problems I have observed that $1$ generally holds, it might be a coincidence. If someone could explain if I am missing something or if either of them(aforementioned cases) is wrong, it would be really helpful.

The zeros of the auxiliary polynomial give you the poles of equal magnitude and opposite signs (either on the real axis or on the imaginary axis). The derivative of the auxiliary polynomial gives you the coefficients that are used in the row with all zero coefficients. It's been a while, but I learned that from reading Modern Control Engineering by K. Ogata (and I had to refresh my memory for answering this question).


How to compute the statistics of the DFT of correlated noise? I would like to know how to compute the statistics of the discrete Fourier transform of a noise signal. To illustrate what I mean, I will first explain in detail a computation I have managed to do myself.

Suppose we have a discrete time series of values $x_n$ with $n$ from 0 to $N-1$. Each $x_n$ is a random variable, uncorrelated with the others, and Gaussian distributed with width $\sigma$. If I define the discrete Fourier transform

$$X_k = \frac{1}{N}\sum_{n=0}^{N-1} x_n e^{-2 \pi i n  k / N}$$

then I find that $X_k$ is a complex random variable with real and imaginary parts Gaussian distributed with width $\sigma/\sqrt{2 N}$. I did the computation by using the fact that the distribution of a sum is the convolution of the distributions, etc.

Now I want to know how to do this computation in the case that $x_n$ are correlated. How does one approach this problem? I can make the assumption that the process is Markovian.

I had originally asked this on the Computation Science site, but I think here is a better fit.

After many weeks I give the answer to my own question.

There is a limit in which we can solve this problem in a reasonably simple way. Suppose we sum enough points in our DFT that the central limit theorem guarantees that the distribution of the sum's real and imaginary parts are Gaussian distributed. Then we only need to compute the variance. If we specialize the case of the variances of the real part of the DFT we can write

$$
\langle (\textrm{Re} X_k) (\textrm{Re} X_l) \rangle = \frac{1}{N^2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}\langle x_n x_m \rangle \cos(2\pi n k / N)\cos(2 \pi m l / N)
$$
The thing to note is that $\langle x_n x_m \rangle$ is just the correlation function of the time domain samples, which we can denote $\rho(n-m)$. Putting this in, we get

$$
\langle (\textrm{Re} X_k) (\textrm{Re} X_l) \rangle = \frac{1}{N^2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1} \rho(n-m) \cos(2\pi n k / N)\cos(2 \pi m l / N)
$$

This form in terms of the correlation function is useful because the correlation function is frequently known; By the Wiener-Khinchin theorem the correlation function is the Fourier transform of the spectral density.

This sum can be computed numerically, or even analytically for some particular forms of $\rho$. To compute covariance with imaginary parts of $X$ just put $\sin$ instead of $\cos$.

I found the idea for this in J. Schoukens and J. Renneboog, IEEE Transactions on Instrumentation and Measurement, Vol. IM-35, No. 3, September (1986).


MATLAB: Adding Noise with Regard to Signal to Noise Ratio (SNR) or EbNo? A signal x takes values from an alphabet set.
My objective is to estimate the channel parameters hest using adaptive methods and then plot the graph of mean square error (MSE) between the estimated coefficients and the actual vs range of SNR values considered. These are the following questions :

Plot of BER vs. signal to noise ratio -- When I want to plot the BER, and the signal takes values other than BPSK say QAM then also I will be estimating the channel coefficients for each snr by adding snr using awgn() using the same code. Then should the X axis be labelled as EsNo or EbNo? In a similar question asked here Adding AWGN noise with a correct noise power to the signal

the answer is to use EbNo. In a Matlab implementation https://www.mathworks.com/matlabcentral/fileexchange/39011-ber-comparison-of-m-ary-qam?focused=5251697&tab=function, using the EbNo range the noise signal of a particular snr is generated. But, the BER plot shown has EbNo on X axis. This is confusing since the awgn() function is using snr and not EbNo.

Plot of error vs snr -- For each noise added, I will calculate the error between the estimates and the actual coefficient. The X axis would be the range of snr and Y Axis the error values. Would the X Axis be snr or EbNo or EsNo irrespective of the symbol set?
In my case the data takes values in symbols and not bits. What should I use? Please help in clearing these concepts. Thank you.

These stuffs $E_s/N_0, E_b/N_0 \textrm{ and SNR}$ are convertible.

\begin{align}
E_s/N_0 &= E_b/N_0 + 10\log_{10}(k) \\
E_s/N_0 &= 10\log_{10}(T_{sym}/T_{samp}) + \mathrm{SNR}
\end{align}

where $k$ is the number of information bits per symbol, $T_{sym}$ is the signal's symbol period and $T_{samp}$ is the signal's sampling period. More details can be found at AWGN channel MATLAB.

And by description of awgn() function:

  y = awgn(x,snr) adds white Gaussian noise to the vector signal x. The
  scalar snr specifies the signal-to-noise ratio per sample, in dB. If x
  is complex, awgn adds complex noise.

In your case, no oversampling thus awgn() uses its snr parameter to generate $E_s/N_0 =$ snr.

You are free to plot your BER curves over $E_s/N_0$ or $E_b/N_0$ because they are equivalent. You just need to convert them back and forth via the modulation order $k$.

However $E_b/N_0$ is usually prefered because (EbN0 wikipedia)

  $E_b/N_0$ directly indicates the power efficiency of the system without
  regard to modulation type, error correction coding or signal bandwidth
  (including any use of spread spectrum). This also avoids any confusion
  as to which of several definitions of "bandwidth" to apply to the
  signal.


Autocorrelation question I am trying to get an understanding of autocorrelation and I am having some issues  with trying to understand the process.

I have a Bernoulli process called $X[t]$. In this process, $P(X[t] = 1) = p$ and $P(X[t] = 0) = 1-p.$

We have a new process formed from $X[t]$ that is:

$$
Y[t] = \frac 12 \left(X[t] + X[t-1]\right)
$$

I have to find the autocorrelation $r_Y [t+\tau,t]$. Here's what I got so far:

\begin{align}
r_Y [t+\tau,t]& = E\left(Y[t+\tau]Y[t]\right)\\
&= \frac 14 E\left\{\left(X[t+\tau] + X[t+\tau-1]\right)\left(X[t] + X[t-1]\right)\right\}\\
&= \frac 14 \left(r_X[t+\tau,t]+r_X[t+\tau,t-1]+r_X[t+\tau-1,t]+r_X[t+\tau-1,t-1]\right)
\end{align}
From there, I can get $r_x[t+\tau,t]$:

$$
r_X[t+\tau,t] =
  \begin{cases}
    p^2       & \quad \text{if } \tau \neq 0\\
    p  & \quad \text{if } \tau = 0\\
  \end{cases}
$$

But the wall that I am running into (and maybe it's because it's not making sense to me because of lack of sleep) is how you substitute that expression back into $r_Y[t+\tau,t]$? I know you end up with three different answers, but trying to substitute it back in seems to cause a mess for me. I am hoping someone will help me with that part.

It is helpful to realize that the auto-correlation of $X[t]$ does not depend on $t$, just on the time difference $\tau$. So we can write it as $r_X(\tau)$. Your result, which is correct, can then be written as

$$E(Y[t+\tau]Y[t])=\frac14\left(r_X(\tau)+r_X(\tau+1)+r_X(\tau-1)+r_X(\tau)\right)\tag{1}$$

Using the expression for $r_X(\tau)$, which is correctly given in your question, $(1)$ can be simplified to

$$E(Y[t+\tau]Y[t])=\begin{cases}\frac{p}{2}(1+p),&\tau=0\\
\frac{p}{4}(1+3p),&\tau\in\{-1,1\}\\p^2,& \text{otherwise}\end{cases}$$

which is of course also independent of $t$.


how to lower frequency of a sound file in matlab? I am new to programming and sound processing,
I have been using matlab to process sound files, a few of which are ultrasonic.
my sampling rate is 250kHz.
I Would like to be able to play the sounds at an audible frequency, there for i need to lower the pitch of the sounds without changing the speed.
The frequency in the sound files are around 70 kHz, and i need to lower them to about 20 kHz (human hearing range 20 Hz - 20 kHz), about 4 times lower.
If i simply change the sampling rate to 250000*4, i get the right frequency but the file become shorter/faster.

I have tried a few things:

Delete every 4th bin of the original vector and multiplying sampling rate by 4. it did the opposite of what i intended and also changed the length.
Duplicate every bin the the audio 4 times and keep the sampling rate, it actually lowered the frequency and the length did not change, but the frequencies were mirrored it 4 times
A method i found online,which i don't understand, this method gave the best result: lower frequency and same file length but it changed the structure of the sound file (see link to spectrogram)

here is the code sample of the last method:

function [signal]=lower_the_freq(original)
% original = rawClip.values';
    binlen=length(original);
    sfq=250000;
    duration=binlen/sfq;

% plot original signal
     subplot(211);

 fourierTransform = fft(original);
 spectrogram(original,1024,512,1024,sfq,'yaxis');    %run the spectrogram
 title('Original signal')

% downsample spectrum by a factor of 2
     n = 2; % downsampling factor
     newSpectrum = fourierTransform(1:n:end);

% zero-pad the positive and negative ends of the spectrum
    pad = floor(length(fourierTransform)/4);
    fourierTransform = [zeros(1,pad) fftshift(newSpectrum) zeros(1,pad)];

% inverse transform
     signal = ifft(length(original)*fftshift(fourierTransform),'symmetric');

% plot the downshifted signal
     subplot(212)
     spectrogram(signal,1024,512,1024,sfq,'yaxis');    %run the spectrogram
     title('Shifted signal')

image

In the end, i took the advice of @chipaudette
in his answer to a similar question
link to answer
and i used Audacity to lower frequency, worked like a charm. in case anyone want's to see the final product, i add a link.
link to manual
link to final product


Signal Plus Weakly Stationary Noise I was reading the book "Spectral Analysis of Time Series" By Herman Koopmans.
On Page 55, he explains that a specific type of non-stationary signal which is the result of adding weakly stationary ergodic noise to a deterministic signal can be decomposed to Wiener spectra. I wonder why there is a need for signal to be ergodic and why weakly stationarity is not enough to derive spectral summation formula for noise and deterministic part?

More specifically based on those assumptions he shows that one has:
$$C_X(\tau) = \lim\limits_{
\tau\to\infty}\frac{1}{2T}\int^{T}_{-T}
X(t + \tau)X(t) dt = C_S(\tau) + C_N(\tau),
\text{ almost surely}.$$

And then
$$F_X(A)=F_S(A)+F_N(A)$$
where $F_Z$ is the spectral distribution of stochastic process $Z(t)$, where
$$X(t)=N(t)+S(t)$$

For ergodic processes, time averages (defined by integrals over time) and ensemble averages (defined by expectations with respect to probability distributions) are identical. This means that the autocovariance is the same, no matter if defined by a time integral or by an expectation:

$$C_X(\tau)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}X(t+\tau)X(t)dt=
E\{X(t+\tau)X(t)\}$$


Fourier series coefficient of signal when Time period is twice the fundamental period My try:

First of all I tried observing the symmetry but I did'nt find any.So I tried to calculate the fourier series coefficient of the signal like this

First I differentiated the signal $x(t)$ so that $y(t)=\frac{dx(t)}{dt}$ Now suppose $d_k$ is the Fourier series coefficient of the differentiated signal and $a_k$ is the fourier series coefficient of Orignal signal $x(t)$
Now $d_k=(jk \omega_0 )c_k$ so $c_k=\frac{d_k}{jk \omega_0}$.Now calculating $d_k$

That is

\begin{align}
d_k &=\frac{1}{3} \int_{-1}^2 \big(\delta(t+1)+\delta(t)-2\delta(t-1)\big) e^{\frac{-jk2\pi t}{T}}dt\\
&=\frac{1}{3}\big(1+e^{\frac{jk2\pi}{3}}-2e^{\frac{-jk2\pi}{3}}\big)
\end{align}

Now $$a_k=\frac{\frac{1}{3}\big(1+e^{\frac{jk2\pi}{3}}-2e^{\frac{-jk2\pi}{3}}\big)}{jk \omega_0}$$
$$a_k=\frac{\big(1+e^{\frac{jk2\pi}{3}}-2e^{\frac{-jk2\pi}{3}} \big)}{j2\pi k}$$

Now checking at

$k=1,a_k\ne 0$
and at $k=2,a_k\ne 0$

so for one odd and one even value of n its not zero so no options matching.What is the mistake?

EDIT:

As stated by @Matt L in the commment to check for period T=6 also, so I did like this

I differentiated the signal from -3 to +3 so the fourier series coefficient I got like this

\begin{align}
d_k &=\frac{1}{6} \int_{-3}^3 \big(\delta(t+3)-2\delta(t+2)+\delta(t+1)+\delta(t)-2\delta(t-1)+\delta(t-2)\big)e^{\frac{-jk2\pi t}{6}}dt\\
&=\frac{1}{6}\big(1+e^{jn\pi}-2e^{\frac{jk2\pi}{3}}+e^{\frac{jn\pi}{3}}-2e^{\frac{-jk\pi}{3}}+e^{\frac{jk2\pi}{3}}\big)
\end{align}

Now at $k=1,d_k=0$ at $k=2,d_k\ne 0$,at $k=3,d_k=0$ but at $k=6$ also $d_k=0$

Now whats the mistake???

As Matt stated, by your calculation you have excluded options a and b but you should also check for options c and d, from which you would see that the answer is the option d.

You could practically check the result of coefficients $d_k$ from the interpretation that CTFS coefficients of a periodic waveform $\tilde{x}(t)$ can be obtained from an inverse-DTFT (or an inverse-DFT more practically), by treating $\tilde{x}(t)$ as a DTFT waveform. Note that the computation is exact if the CTFS coefficients $d_k$ were of finite length in advance, but otherwise approximate which would only get better by choosing long enough samples to represent the waveform $\tilde{x}(t)$.

Therefore, the following computes (approximately here) $3N$ CTFS coefficients of $x(t)$:

>     N = 128;
>     x = [ones(1,N) , -1*ones(1,N) , zeros(1,N)];
>     d = ifft(x);

On the orther hand the following code would do the same by treating $\tilde{x}(t)$ as periodic in $2 T_0$ rather than $T_0$

>     N = 128;
>     x = [ones(1,N) , -1*ones(1,N) , zeros(1,N)];
>     d = ifft([x x]);

Based on your edit, the following provides you the answer:
Given that the CTFS coefficients of the derivative signal is:

\begin{align}
d_n &=\frac{1}{6} \int_{-3}^3 [\delta(t+3)-2\delta(t+2)+\delta(t+1)+\delta(t)-2\delta(t-1)+\delta(t-2)]e^{\frac{-jn2\pi t}{6}}dt\\
&=\frac{1}{6}[  e^{j\frac{2\pi}{6}3n} - 2 e^{j\frac{2\pi}{6}2n} + e^{j\frac{2\pi}{6}1n} + 1 -2e^{-j\frac{2\pi}{6}1n} + e^{-j\frac{2\pi}{6}2n}] \\
&=\frac{1}{6} [  e^{j\pi n} - 2 e^{j\frac{2\pi}{3}n} + e^{j\frac{\pi}{3}n} + 1 -2e^{-j\frac{\pi}{3}n} + e^{-j\frac{2\pi}{3}n}] \\
\end{align}

Above three lines were straightforward. Now we shall group those terms to yield something that can be simplified:

$$d_n =\frac{1}{6} [  (1 + e^{j\pi n}) - 2 ( e^{j\frac{2\pi}{3}n} + e^{-j\frac{\pi}{3}n} ) + (e^{j\frac{\pi}{3}n} + e^{-j\frac{2\pi}{3}n})] $$

Now add $2\pi n$ to the negative angles so that they become: $$e^{-j\frac{\pi}{3}n} = e^{j( 2\pi n -\frac{\pi}{3}n) } = e^{j\frac{5\pi}{3}n } = e^{j\pi n }e^{j\frac{2\pi}{3}n } $$ and  $$e^{-j\frac{2\pi}{3}n} = e^{j( 2\pi n -\frac{2\pi}{3}n) } = e^{j\frac{4\pi}{3}n } = e^{j\pi n}e^{j\frac{\pi}{3}n } $$ respectively. And plugging them into $d_n$ line, yields:

\begin{align}
d_n &=\frac{1}{6} [  (1 + e^{j\pi n}) - 2 ( e^{j\frac{2\pi}{3}n} + e^{j\pi n }e^{j\frac{2\pi}{3}n } ) + (e^{j\frac{\pi}{3}n} + e^{j\pi n}e^{j\frac{\pi}{3}n})]\\
&=\frac{1}{6} [  (1 + e^{j\pi n}) - 2 (1 + e^{j\pi n}) e^{j\frac{2\pi}{3}n} + (1 + e^{j\pi n}) e^{j\frac{\pi}{3}n}]\\
&=\frac{1}{6} [1 + e^{j\pi n}]\cdot[1 - 2 e^{j\frac{2\pi}{3}n} + e^{j\frac{\pi}{3}n}]\\
&=\frac{1}{6} [e_n]\cdot[f_n]\\
\end{align}

Now it can be shown that the product term $e_n$ will be zero for all odd indice $n$, hence $d_n$ will also be zero for all odd $n$. Note that the other product term $f_k$ will be zero for $n = 6m$ that $d_n$ will also be zero for $n = 6,12,18...$ However for your example it sufficies to show that $e_n$ will be zero for odd $n$.

FURTHERMORE below is the theoretical proof that it's not a coincidence to have all the odd indexed terms in CTFS to be zer0, independent of the signal itself, when the period $T_y$ of the periodic signal $x(t)$ is assumed to be twice that of its fundamental period $T_x$. Assume $x(t)$ is a periodic signal with fundamental period of $T_x$ and let $y(t)$ be the same signal $x(t)$ but interpreted to have a period of $T_y = 2 T_x$. Let $a_k$  and $b_k$ denoted the CTFS coefficients of $x(t)$ and $y(t)$ respectively, then we have:

\begin{align}
b_k &= \frac{1}{T_y} \int_{0}^{T_y} y(t) e^{-j\frac{2\pi}{T_y} k t} dt = \frac{0.5}{T_x} \int_{0}^{2T_x} x(t) e^{-j\frac{2\pi}{T_x} (k/2) t} dt\\
&= 0.5 \left( \frac{1}{T_x} \int_{0}^{T_x} x(t) e^{-j\frac{2\pi}{T_x} (k/2) t} dt + \frac{1}{T_x} \int_{T_x}^{2T_x} x(t) e^{-j\frac{2\pi}{T_x} (k/2) t} dt \right)\\
\end{align}

Make the substitution $t' = t - T_x$ in the second integral and recognize that $x(t'+T_x) = x(t')$ as $x(t)$ is periodic with $T_x$ , yielding:

$$b_k = 0.5 \left( \frac{1}{T_x} \int_{0}^{T_x} x(t) e^{-j\frac{2\pi}{T_x} (k/2) t} dt + e^{-j\frac{2\pi}{T_x} (k/2) T_x} \frac{1}{T_x} \int_{0}^{T_x} x(t') e^{-j\frac{2\pi}{T_x} (k/2) t'} dt' \right) $$

Now recognize the integrals as $a_{k/2}$; i.e, the CTFS of the signal $x(t)$ evaluated at $k/2$. And simplify the sum:
\begin{align}
b_k &= 0.5 \left( a_{k/2}+ e^{j\pi k} a_{k/2} \right)\\
&= 0.5 \left( 1+ e^{j\pi k} \right) a_{k/2} \\
&= \begin{cases}{ a_{k/2} ~~, \text{ for k even, k=2m  , m=1,2,...}\\ 0 ~~~~~~~, \text{ for k odd  , k=2m+1, m=1,2,...}}\end{cases} \\
\end{align}


How Can JPEG Be Used to Realize a Transfer Application with a Fixed Bit Rate? For simplicity, assume a fixed bandwidth (5MB/s) and that each second an image with 12MP is generated. The quality of the image should be as high as possible (in terms of PSNR). Which compression rate is necessary and how would the task differ or become simplified when using JPEG 2000 instead?

Thoughts

When the image is a colored image with 24bpp, then the size of a 12MP image would be 12*10^6*24/8 Byte = 36 MB. That divided by 5MB yields a compression rate of approximately 7. Then transferring works as desired. However, I don't see how the task changes when using JPEG 2000.

You calculated the needed compression ratio correctly.
Now to answer the rest of the question:

The JPEG Algorithm doesn't let you set the size of the output file only the Quality measure. Namely to ensure the size of the output is close to the limit (Namely highest quality) requires trial and error (Though using the Quality you can get a measure of the bits per color).
The JPEG 2000 on the other hand allows you to cut the data at the exact number of bytes required hence makes the task much simpler (You can practically set a threshold on the output size).
The image quality of JPEG at compression ratio of ~100 is really bad (See Wikipedia article on JPEG). While JPEG 2000 can get better results.

To understand the Quality in JPEG, have a look at How JPG Works by Colt McAnlis:

Pay attention to the line:

  Probably most important, is that the quality parameter varies depending on the image. Since each image is unique, and presents different types of visual artifacts, the Q value will be unique as well.

This is what JPEG 2000 makes easier.


Subtracting signals to re-construct the middle Below image represents a sample of some of my captured data. There are three waves if you will. The middle wave has a very large dip ( cluttered by noise ) .

My goal here is to try and resurrect the middle wave to be noise free ( as much as possible ) . I will know I succeeded when the middle wave looks very 'similar' to the other two .

Here is what I am thinking, but I am not sure if the idea I have in mind is mathematically correct.

Given the red portion ( samples between 1920 to 2080 ) and the yellow portion ( samples between 2020 to 2170 ), I can take two FFT's. I can then determine the frequency that is dominant in both sample sets. If I take the average between these two ( FFT from red portion and FFT of yellow portion ), would that give me a mathematically correct estimation of what the middle wave should be?

No, you can't restore a low - frequency signal by DFT of its parts due to the periodic DFT nature and because
$$
IDFT(DFT( \frac{U(t) + V(t)}{2})) == \frac{U(t) + V(t)}{2}
$$
by definition.
You can apply filtering or smoothing to remove the high - frequency noise.


Given the absolute value of a function, find hilbert pairs that make it up This might be a solved problem (or a known unsolved problem) because it is related with finding a causal signal that has a desired magnitude response.

$$F(f)=R(f)+iI(f)$$

$$F(f)=R(f)+i \int^\infty_{-\infty}\frac{R(u)}{f-u}du$$

$$F(f)F^*(f)=R^2(f)+\bigg(\int^\infty_{-\infty}\frac{R(u)}{f-u}du\bigg)^2$$

We only know $F(f)F^*(f)$

Is it possible to solve this, is the solution unique?

This problem is known as spectral factorization. For convenience, let me use $j\omega$ as the function argument. Assuming that $F(j\omega)$ is rational, i.e. realizable by lumped elements (R,C,L), the squared magnitude $F(j\omega)F^*(j\omega)$ is a non-negative real rational function. In the Laplace transform domain, this function is given by

$$G(s)=F(s)F^*(-s^*)\tag{1}$$

because $F^*(j\omega)$ corresponds to the time domain function $f^*(-t)$, the Laplace transform of which is $F^*(-s^*)$. From $(1)$ it is clear that for every pole $s_{\infty}$ of $G(s)$ there must be a pole at $-s_{\infty}^*$, i.e. mirrored across the imaginary axis. The same is of course true for the zeros.

Now you can factor $G(s)$ in such a way that you assign all poles to the left of the imaginary axis (i.e. with negative real parts) to $F(s)$. If causality is your only requirement, you are free in assigning the zeros, apart from the need to always take complex conjugate pairs, which guarantees that the impulse response $f(t)$ corresponding to $F(s)$ is a real-valued function. By assigning all poles with negative real parts to $F(s)$, you make sure that $F(s)$ is a causal and stable system. If $G(s)$ has poles on the imaginary axis, there exists no stable and causal system $F(s)$.

So the answer to your question is: yes, it is a known problem and it can be solved; no, the solution is not unique (because you can freely choose which zeros you assign to $F(s)$).


Zero, First, Second ... nth-order Hold The rectangular function is defined as:
$$\mathrm{rect}(t) = \begin{cases}
0 & \mbox{if } |t| > \frac{1}{2} \\
\frac{1}{2} & \mbox{if } |t| = \frac{1}{2} \\
1 & \mbox{if } |t| < \frac{1}{2}. \\
\end{cases}$$

The triangular function is defined as:
$$\operatorname{tri}(t)
=
\begin{cases}
1 - |t|, & |t| < 1 \\
0, & \mbox{otherwise}
\end{cases}
$$
It is the convolution of two identical unit rectangular functions:
$$
\operatorname{tri}(t) = \operatorname{rect}(t) * \operatorname{rect}(t) \quad
=  \int_{-\infty}^\infty \mathrm{rect}(\tau) \cdot \mathrm{rect}(t-\tau)\ d\tau$$
Zero-order hold and First-order hold use these functions. In fact, it has:
$$
x_{\mathrm{ZOH}}(t)\,= \sum_{n=-\infty}^{\infty} x(n)\, \mathrm{rect} \left(t-n\right) \
$$
for Zero-order hold, and
$$x_{\mathrm{FOH}}(t)\,= \sum_{n=-\infty}^{\infty} x(n)\, \mathrm{tri} \left(t-n\right) \ $$
for First-order hold. Since $\operatorname{tri}(t) = \operatorname{rect}(t) * \operatorname{rect}(t)$, I would like to know if this is just a coincidence or if, for the Second-order hold the impulse response is
$$\operatorname{tri}(t)*\operatorname{tri}(t)=\left(\operatorname{rect}(t) * \operatorname{rect}(t)\right)*\left(\operatorname{rect}(t) * \operatorname{rect}(t)\right).$$
Is it true also for a general $k$-th order hold? Namely, put
$$x_{\mathrm{K-TH}}(t)\,= \sum_{n=-\infty}^{\infty} x(n)\, \mathrm{g}_k \left(t-n\right) \ $$
where $\mathrm{g}_k \left(t-n\right)$ is the impulse response of the $k$-th order hold, I would like to know if its impulse response is
$$\mathrm{g}_k \left(t-n\right)=\left(\operatorname{rect}(t) * \operatorname{rect}(t)\right)* \dots * \left(\operatorname{rect}(t) * \operatorname{rect}(t)\right),$$
k times.

This is not the case. First of all, a second-order hold would use three sample points to compute an interpolation polynomial, but your suggested impulse response $\text{tri}(t)\star\text{tri}(t)$ is non-zero in an interval of size $4$ (assuming a sample interval of $T=1$, as you do in your question). However, the impulse response corresponding to a second-order hold must have a support of length $3$.

Now you could suggest that an $n^{th}$-order hold could have an impulse response which is the convolution of $n$ rectangular functions. In this case you would get the correct support size, but of course that's not sufficient.

An $n^{th}$-order hold computes a piece-wise interpolation using $n+1$ consecutive data points. This is in analogy with a zero-order hold using a single data point, and a first-order hold, which uses two data points. This definition is commonly used in the literature (see e.g. here and here).

It's straightforward to show that the second-order polynomial that interpolates three data points $y[-1]$, $y[0]$, and $y[1]$ is given by

$$P(t)=y[-1]\frac{t(t-1)}{2}+y[0](1-t^2)+y[1]\frac{t(t+1)}{2}\tag{1}$$

In order to find the impulse response achieving the interpolation given by $(1)$, we have to equate $(1)$ with the expression

$$y[-1]h(t+1)+y[0]h(t)+y[1]h(t-1)\tag{2}$$

If we choose the support of the impulse response $h(t)$ as the interval $[-1,2]$, which is equivalent to choosing the interpolation interval $[0,1]$, equating $(1)$ and $(2)$ results in the following impulse response of a second-order hold:

$$h(t)=\begin{cases}\frac12(t+1)(t+2),&\quad -1<t<0\\
1-t^2,&\quad0\le t \le 1\\
\frac12 (t-1)(t-2),&\quad 1<t<2\\
0,&\quad\text{otherwise}\end{cases}\tag{3}$$

The impulse response $(3)$ of a second-order hold looks like this:

I leave it up to you to show that this impulse response cannot be generated by convolving three rectangular functions with each other.


scanning an image for projected points I'm projecting a series of points onto a surface and then recording their locations in the corresponding image.  However, if one goes pixel by pixel a single point may, at times, take up multiple pixels and would thus be recorded multiple times.  What is the best (least error prone) way to determine the location of a point?  I've thought of simply scanning the image line by line checking to see if a pixel corresponds to that of a point, grouping the pixel clusters together based upon distance, and taking the center of this cluster to be the coordinate of the point.  Perhaps something like an R-tree where I can order things by distance?

Any suggestions? Or tried and true techniques?

Assuming that they are adjacent and have a "red" energy/pixel value distinct from not projected areas, you may use a simple 2D filter, that fits the size of expected dispersion, and that way get the local maximum that represents the center of the collection of points that represents a single projection.


Need an Identity for CTFT Polynomial Raised to an Exponent I need to the find the inverse continuous time Fourier transform for unitary angular frequency of the following signal:

$e^{a\omega^2 - b\omega + c}$ where $a$ and $b$ and $c$ are real numbers and I want to find the corresponding time-domain signal.

I was able to get the $e^{a\omega^2 + c}$ term transformed, but I was unable to find a transform for $e^{- b\omega}$ that didn't result in a distribution.

Without doing math and looking at the symmetry, I want to say that the $e^{- b\omega}$ term does a complex time shift, but now my mind hurts. Can somebody look up the identity?

--edit--
This spectrum doesn't need to be symmetric, or causal, so expecting some solutions to have imaginary time could be correct. Can somebody comment on this?

Dilip's idea is an excellent suggestion. Let me slightly modify it. You can rewrite
$$a\omega ^{2}-b\omega +c = a(\omega -\omega_{0} )^{2}+d$$
So the whole thing can be written as

$$e^{d}\cdot e^{a(\omega -\omega_{0} )^{2}}$$

That's essentially a Gaussian shifted in frequency with a constant scale factor. The transform of a Gaussian is just another Gaussian and the shift in one domain corresponds to multiplication with a complex exponential in the other domain so the inverse ought to look something like:

$$x(t) = k\cdot e^{\tilde{a}t^{2}}\cdot e^{j\omega _{0}t}$$
where $k$, $\tilde{a}$, and $\omega _{0}$ are determined by the original a, b, and c. a better be a negative number, though, otherwise I don't think that the transform would be convergent.


What is the purpose of the filtering the image data in the PNG file format I am reading the W3C description of the PNG file format, in particular section $9$ on filtering. The filter functions:

None(x) : the identity filter
Sub(x) : subtract the pixel value of $x$ by the pixel value in the position immediately before x
Up(x):   subtract the pixel value of $x$ by the pixel value in the scanline before the scanline containing $x$
Average(x): Subtract the pixel value of $x$ by the mean of the subtracted values described in Su and Up
Path(x) : Substract the pixel value of $x$ by the result of the Paeth Predicator at three values, see the spec.

are applied to different scanlines of an image; the filter function is applied to each pixel in a line.
It is stated without justification that:

Filtering transforms the PNG image with the goal of improving compression

My question is why does applying these filters improve compression?

The usage of a word filtering here is somewhat of a misnomer. The filtering operation is typically associated with waveform processing for various applications, as pulse shaping in the transmitter or noise suppression in the receiver. Image filters are also used, but primarily not for better compressibility of image files. In the context of your reference of PNG file format standard, it is possible that this operation would be more aptly named data preprocessing.
A Wikipedia article on PNG file format is more explicit on the compressibility of filtered PNG image data:

An image line filtered in this way is often more compressible than the
raw image line would be, especially if it is similar to the line
above, since the differences from prediction will generally be
clustered around 0, rather than spread over all possible image values.

To further clarify the effect of PNG data filtering on data compressibility, I cite the LZ77 algorithm article:

LZ77 algorithms achieve compression by replacing repeated occurrences
of data with references to a single copy of that data existing earlier
in the uncompressed data stream.

The less variability of "the differences from prediction" that are "clustered around 0", the more "repeated occurrences of data" exist "in the uncompressed data stream". Therefore, the more effective is the LZ77 algorithm's work.
The LZ77 is an algorithm the PNG standard uses for compression.
Notice also the statement from  W3C description of the PNG file format:

Filters are applied to bytes, not to pixels, regardless of the bit
depth or colour type of the image.

Again, the Wikipedia article is helpful to better understand this statement w.r.t. compressibility:

An image line filtered in this way is often more compressible than the
raw image line would be, especially if it is similar to the line
above, since the differences from prediction will generally be
clustered around 0, rather than spread over all possible image values.
This is particularly important in relating separate rows, since
DEFLATE has no understanding that an image is a 2D entity, and instead
just sees the image data as a stream of bytes.


Segmenting an image into columns I want to segment the image below (click for bigger image) into columns based on the natural gaps that exist in the columns. The columns should be perfectly vertical and need not be the same size, but should split each column by itself.

I have tried Hough Line Transform, Scipy Label, convolution filters, specifying manual cuts, etc. all with very limited success and nothing I would consider a "solution".

What approach should I take to segment these types of images into columns? Preferably a solution that I can code myself or exists in OpenCV or Python.

You may not need to delve into scipy or opencv.  The following solution uses pil and numpy:

from PIL import Image
import numpy

#define a greyscale threshold which we will use later
threshold = 230

#open image as greyscale with PIL and convert to numpy array
with Image.open('your_image.png').convert('L') as img:
    arr_gry = numpy.array(img)

#create 1-D array of max column values
colmax = numpy.zeros(arr_gry.shape[0])
colmax = numpy.min(arr, axis=0)

#analyse the bands array and define start and end of each band
band_limits = {}
band_list = []
prev_column = 0
for i,column in enumerate(colmax):
    if column > threshold and ('start' not in band_limits):
        band_limits['start'] = i
    if column <=threshold and prev_column > threshold:
        band_limits['end'] = i
        band_list.append(band_limits)
        band_limits = {}
    prev_column = column

###The following is for illustration purposes only

#convert to RGB array
arr_clr = numpy.stack([arr_gry,arr_gry,arr_gry], axis=2)

#Colour the bands and mid points
for band in band_list:
    arr_clr[:,band['start']:band['end']] = [255,0,0]
    arr_clr[:,(band['start']+band['end'])//2] = [0,0,255]

#Convert array to Image
out = Image.fromarray(arr_clr)
out.show()

This produces the following image

Obviously you can use opencv if desired to open, view and save the images.  But in this example all the processing is done in numpy so pil is more than capable.


How do I implement the Structure tensor (in Matlab)? I am trying to figure out the details on how to implement the structure tensor in Matlab and need some advice!

For an image $\ I(x,y) $ the structure tensor S is given by:

$$ S=\begin{pmatrix}
W \ast I_x^2 & W \ast (I_xI_y)\\
W \ast (I_xI_y) & W \ast I_y^2
\end{pmatrix}$$
where W is a smoothing kernel, $\ \ast $ denotes convolution and subscript denotes partial derivative with respect to.

I want to use a Gaussian as smoothing kernel. This is separable:
$$ W = G_{\sigma}*(G_{\sigma})^T$$
The constant $\ \sigma $ controls the outer/integration scale.
E.g. if I truncate it at $\ +/- 2\sigma$ and use  $\ \sigma = 1$ I get:

$$ G_{\sigma} = [0.0103,0.2076,0.5642,0.2076,0.0103] $$

However the input to the Gaussian might have a mean different from 0.

How do I compensate for this?

I was thinking substracting the mean first.

I could use central differences to estimate the partial derivatives:
$$ I_x \approx [-0.5 \hspace{2 pt} 0 \hspace{2 pt} 0.5]*I $$
However since the image is noisy it is probably better to use a more robust differentian operator.
The Sobel operator is such an operator. However according to wikipedia it is not perfectly rotational symmetric. The Scharr operators are supposed to be perfectly rotational symmetric and robust approximations of the partial derivatives:

$$ I_x \approx \frac{1}{32}\begin{pmatrix}
-3 & 0 & 3\\
-10 & 0 & 10\\
-3 & 0 & 3
\end{pmatrix} * I$$

$$ I_y \approx \frac{1}{32}\begin{pmatrix}
-3 & -10 & -3\\
0 & 0 & 0\\
3 & 10 & 3
\end{pmatrix} * I$$

The Scharr operators are separable:
$$ \begin{pmatrix}
-3 & 0 & 3\\
-10 & 0 & 10\\
-3 & 0 & 3
\end{pmatrix} = [-1,0,1]*[3,10,3]^T $$

$$ \begin{pmatrix}
-3 & -10 & -3\\
0 & 0 & 0\\
3 & 10 & 3
\end{pmatrix} = [3,10,3]*[-1,0,1]^T$$

The elements of the structure tensor now becomes series of consecutive horizontal and vertical convolutions with the original image.

EDIT: Another option (thanks to nikie!) is the Gaussian derivative.
Introducing
$$ GD_{\sigma}=\frac{-x}{\sigma}G{\sigma}$$

Is this then correct?
$$ Ix = G{\sigma}^T*(GD_{\sigma}*I) $$
$$ Iy = G{\sigma}*(GD_{\sigma}^T*I) $$

How do I Implement this easily in Matlab?

How about in C++?

Any comments regarding choice of smoothing kernel and differentiation operator?

EDIT 2: My Matlab code so far:

% I: image
% si: inner (differetiation) scale,
% so: outer (integration) scale
function [Sxx, Sxy, Syy] = structureTensor(I,si,so)
[m n] = size(I);

Sxx = NaN(m,n);
Sxy = NaN(m,n);
Syy = NaN(m,n);

% Robust differentiation by convolution with derivative of Gaussian:
x  = -2*si:2*si;
g  = exp(-0.5*(x/si).^2);
g  = g/sum(g);
gd = -x.*g/si; % is this normalized?

Ix = conv2( g',conv2(gd,I) );
Iy = conv2( g,conv2(gd',I) );

Ixx = Ix.^2;
Ixy = Ix.*Iy;
Iyy = Iy.^2;

% Smoothing:
x  = -2*so:2*so;
g  = exp(-0.5*(x/so).^2);
Sxx = conv2( g',conv2(g,Ixx) );
Sxy = conv2( g',conv2(g,Ixy) );
Syy = conv2( g',conv2(g,Iyy) );

How do I normalize the Gaussian derivative?

However the input to the Gaussian might have a mean different from 0.

Unless your image is completely constant, $I_x^2$ and $I_y^2$ are guaranteed to have a mean > 0 (assuming the input pixels are all reals, of course).

  How do I compensate for this?

Why would you want to do that?

  I was thinking substracting the mean first.

Subtracting the mean before applying a derivative operator like Sobel or Scharr? That shouldn't make any difference.

  The Scharr operators are supposed to be perfectly rotational symmetric

Perfect rotational symmetry? Is that even possible with a finite kernel? I think the Scharr operator is just optimal (in a least-squares sense) for that filter size. With larger filter sizes, you get better rotational symmetry.

  How do I Implement this easily in Matlab?

You can use imfilter to apply the Sobel or Scharr filters, then .* to multiply/square the results pixelwise, then imfilter again to apply the gaussian smoothing kernel.

  How about in C++?

Any decent image processing library should contain functions for linear filters and pixel-wise multiplication. e.g. for the IPP, the functions would be called ippiFilterRow/Column..., ippiMul... and ippiSqr...

  Any comments regarding choice of smoothing kernel and differentiation operator?

I think the gaussian smoothing kernel is pretty standard. Sobel and Scharr are both common choices for the differentiation operator. Another good choice would be a gaussian derivative filter. Which one's best depends on your application: how noisy are your images? How sharp are the contrasts you're looking for? (If there's lots of noise, and the contrasts are sharp, a derivative filter with a larger aperture like the gaussian derivative might be better). How important is rotational symmetry?

  How do I normalize the Gaussian derivative?

The gaussian is: $\frac{1}{2 \pi  \sigma ^2}e^{-\frac{x^2+y^2}{2 \sigma ^2}}$

Deriving by x yields:

$\partial _x\left(\frac{1}{2 \pi  \sigma ^2}e^{-\frac{x^2+y^2}{2 \sigma ^2}}\right) = -\frac{e^{-\frac{x^2+y^2}{2 \sigma ^2}} x}{2 \pi  \sigma ^4}$

So I think you should divide by gd = -x.*g/(si^2); because sum(g)=2*pi*(si^2)


Why there is no paper optimize the consumption power of the circuit in energy harvest? In the optimization of  energy harvest ,some paper want to maximize the harvested power ,rate or minimize the transmitted power,so when they write the formula of optimization,they may assume the beamforming,transmitted power  power splitter ,power allocation as variable.

However,i am confused that why there is no one assume the "consumption of circuit" as variable ? All of paper are either assume it is a constant value or just ignore them.

Does anyone know or have a good explanation about we don't assume  "consumption of circuit" as variable?

In my opinion,i think the purpose of optimizing something in the energy harvesting is that we want to let the life of battery become longer when we "transmit" the signal correctly,so we will focus on something which has a relation with the transmitted power,like SNR  or power splitter,beamforming direction,etc.

Does anyone know or have a better explanation about it?Because we do need a circuit to transmit the signal,and the signal does consumpt the power

Everyone wants to minimize the power of their circuits, so there's nothing special about energy harvesting devices, at all.

Also, every paper about making communications more reliable while is about making the power usage lower, simply because if you want constant reliability, you can then achieve that with lower power.

Basically, all the research on transmitters  uses the energy per bit as the critical measure, so you really don't need to look for energy harvesting-specific papers.

So, this really boils down to basic textbook knowledge of transmitter design, not to speciality papers.

Now, that's on the theoretical side of things. When it comes to building silicon doing the things that you derived are good in theory, then you'll find that modern silicon's manufacturing techniques are very much limited by the power consumption of devices – so, the secrets that you're asking for a worth billions of Qualcomm's, Apple's, ZTE's, and (formerly) intel's money.


How do I sum the discrete cosine transform cosine waves? Sorry if this sounds like a stupid question, I'm new to DCT.

I've implemented various types of DCT (like IV and II) and they all give vastly different coefficients. From what I've read, the DCT expresses a signal as a sum of cosine waves with various oscillating frequencies. How can different formulas that give different coefficients express the same signal? I created 7 cosine waves (With one DC coefficient) on a graph with proportional periods according to N, and summed them but end up with a function that doesn't represent the signal at all. I'm probably doing this wrong, but here's what I did:

Calculated the 1D-DCT-II of this array:
 [1 2 3 4 5 6 7 8]
DCT:

[12.727922061358 -6.4423230227051 -1.3322676295502E-15 -0.67345480090394 0 -0.200902903736 -1.5765166949677E-14 -0.050702322759621]

Now I don't know if this is how to do it but I tried this (with a window of length 8):

12.72 + -6.44cos((pi/7)x)+ 0cos((pi/3.5)x)...

I can use the inverse DCT (DCT-III) and get the input values back, but I want to see the signal that represents the original input values. I tried summing the frequencies but I don't know if I'm doing it right so can someone guide me?
Btw the formula I used was

$y[$n]= $realarray[$n]*cos((M_PI/$N)*($n+0.5)*($k));

You used the DCT-II with a certain scaling convention:

$$X[k]=s[k]\sum_{n=0}^{N-1}x[n]\cos\left(\frac{\pi}{2N} (2n+1)k\right),\qquad k=0,1,\ldots,N-1\tag{1}$$

where $s[0]=1/\sqrt{N}$, and $s[k]=\sqrt{\frac{2}{N}}$ for $k=1,\ldots,N-1$. Its inverse is given by an appropriately scaled DTC-III:

$$x[n]=\sum_{k=0}^{N-1}s[k]X[k]\cos\left(\frac{\pi}{2N} (2n+1)k\right),\qquad n=0,1,\ldots,N-1\tag{2}$$

where $s[k]$ are the same scaling coefficients as used in $(1)$.

So, given the scaled DCT-II coefficients $X[k]$, the original signal can be written as

$$x[n]=\sqrt{\frac{2}{N}}\left[\frac{1}{\sqrt{2}}X[0]+X[1]\cos\left(\frac{\pi(2n+1)}{2N}\right)+X[2]\cos\left(\frac{\pi(2n+1)}{N}\right)+\ldots\right]$$


Design narrow bandpass filter for signal with high sampling rate I have to bandpass-filter a signal which has been sampled with 4000 Hz.

Only frequencies around 15 Hz shall remain after filtering (let's say in a band between 10 Hz and 20 Hz; the narrower the better).

I have several questions here:

What is the recommended way to perform tasks like this?
What bandpass filter is suitable?
Do I need to downsample the signal before filtering. If yes, what type of additional low pass filter should I use to avoid aliasing?

In my case, it is important that the phase of the original signal will not be distorted or shifted in the whole process. Furthermore, the filter step response should have close to zero overshoot. The settling time could therefore be a little larger.

I do not want to perform a sophisticated filter design which delivers optimal results for exactly this specific signal. I am more interested in a general solution that delivers solid results using standard IIR or FIR filters (e. g. as available in Python's SciPy library) which could be reused afterwards for similar tasks.

UPDATE:

After the answer of @MarcusMüller and the provided links, I basically tried every FIR filter design method available in Python's SciPy library and deciced to go with remez. I developed the following code which may be used for further discussion:

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import lfilter, remez

F_test = 20.0
duration = 10.0
fs = 8000
samples = int(fs*duration)
t = np.arange(samples) / fs
signal_test = (5.0 * t * np.sin(2.0*np.pi*F_test*t)) + (0.5 * np.sin(2.0*np.pi*5.0*t)) + (0.5 * np.sin(2.0*np.pi*100.0*t))

#design filter

ntaps = 5000
edges = [0, F_test - 5.0, F_test - 2.5, F_test + 2.5, F_test + 5.0, 0.5 * fs]
taps = remez(ntaps, edges, [0, 1, 0], Hz=fs, maxiter=2500)

#apply filter
signal_test_filtered = lfilter(taps, 1, signal_test)

#create plot
fig = plt.figure()
ax0 = fig.add_subplot(111)
ax0.plot(t, signal_test_filtered, label='signal_test_filtered')
ax0.set_xlabel("time [s]")
ax0.legend()
fig.show()

So, first, to put things into perspective: 4kHz is not a high sampling rate these days (add 5 orders of magnitudes, and things become hard).

Your 15 kHz passband doesn't say anything about the complexity of the filter; what counts is the transition width, ie. the distance between pass- and stopband, as well as the attenuation of the stopband. For a slightly specialized answer, see the answers to http://dsp.stackexchange.com/questions/31066/how-many-taps-does-an-fir-filter-need/31077.

  What is the recommended way to perform tasks like this?

Do a filter design, apply filter.

  What bandpass filter is suitable?

Any bandpass filter that suits your needs – which you haven't specified.

What's missing is:

transition width
stop-band attenuation
acceptable ripple
length constraints

What we know is:

You want a linear phase filter (because linear phase means constant group delay, and that means un-broken phase relationships)
Thus, you probably want a FIR that is symmetric in time

  Furthermore, the filter step response should have close to zero overshoot.

Gibb's phenomenon is non-negotiable :) so yeah, use a longish filter with a nice rolloff.

  The settling time could therefore be a little larger.

That will be the effect, yes. You could use a filter design tool that allows you to design with a windowing method and use a window that suits your application well – which I don't know.

  I am more interested in a general solution that delivers solid results using standard IIR or FIR filters (e. g. as available in Python's SciPy library)

Well, yeah, that's what I'd recommend, but you say just shortly above that you don't want to do a filter design? I'm a bit conflicted. Anyway, there's a lot of functions that will give you a proper design.

I'd recommend the following:

Design a low-pass FIR filter
convert it to a bandpass filter, by multiplying with a cosine (real-valued data, symmetric filter) or $e^{j2\pi\frac{f_{center}}{f_s}\cdot n}$ (complex-valued data, one-sided filter) as needed.

Designing is easy, something along (untested, straight from the back of my head)

from scipy.signal import fir_filter_design
from math import cos, pi

f_center=ToBeDefined!!!
f_cut = 15
f_s = 4e3
f_rel = f_s/2/f_cut
transition_width = ToBeDefined!!! # e.g. 10
trans_rel = f_s/2/transition_width
attenuation = ToBeDefined!!! # e.g. 10.0**-5

num_of_taps = estimate_by_reading_my_answer_above(f_rel, trans_rel, attenuation)

lowpass_taps = fir_filter_design.firwin2(numtaps=num_of_taps, cutoff=f_rel, window = "hamming"|"blackmanharris"|"hann"|"chebwin"|…)

bandpass_taps = [lp_tap * cos(2*pi*f_center/f_s*n) for n, lptap in enumerate(lowpass_taps)]

  Do I need to downsample the signal before filtering. If yes, what type of additional low pass filter should I use to avoid aliasing?

Well, considering you won't need most of your signal, yes, that seems advantageous. Typically, you'd just filter away (and decimate in the same step) an 1/N of your initial sampling rate. For example, if your bandpass lies about 100 Hz - 15 / 2 Hz to 100 Hz + 15 / 2 Hz, you won't need anything above 200 Hz – so just decimate to 1/20 of your input sampling rate, with a low pass filter that has 1/40 input sampling rate transition width.

But again, 4kHz is laughably little, and a couple thousand taps won't stress your PC at all. If in doubt,

from scipy import signal
signal.convolve(your_input_signal, filter_taps, method="fft")

does a fast convolution. (you can also just omit the method, or set it to auto, because scipy will just pick fast convolution for long filters, automatically). For a bit of an impression on what's possible with halfway-optimized code on a PC, see my answer to how to implement signal generation on a GPU (and why not); there, I showed that you can do a 107-tap filter at 160 MS/s on my PC back then, so you should be able to do a let's say 700 tap filter at easily 20 MS/s (just a wild guess). That's only 5 thousand times as fast as your sampling.


Implementing a 1-D Kalman Filter Regression, Missing the smoothing action (getting the opposite) I am implementing the 1D Kalman Filter in Python on a fundamentally noisy set of measurement data, and I should be observing a large amount of smoothing...but, instead, my Kalman Filter is doing the exact opposite of smoothing, so that my filtered results look like this:

And since my "model" is just the simple trend formula, $$x_{k} = 3*(x_{k-1} - x_{k-2}) + x_{k-3} $$
the model comes up with outrageous predictions when the expected smoothing doesn't happen:

I am working off an incompletely specified 1D Kalman Filter implementation, so I am having to fill in the gaps...but I should be seeing much softer filtered values, with the peaks and valleys being outliers, with my naive model.  Also, my predictions should actually be more modest than the actual data (which would follow from a smoothed filtered collection).

I am trying to figure out what I am missing.  Here is the pseudo code for what I have implemented so far:

 Step 1: Minimize T = -log(Q/R) in order to discover Q
         Here I basically run the kalman filter with given T values,
         calculating Q from the Residual std

 Step 2: Run Kalman Filter with T held constant

 xprior = [x1,x2,x3] (calculated from the first three values of data)

 #prediction
 x_{k|k-1} = model_estimate(xprior)

 #add measurement
 innovation = data[k] - x_{k|k-1}
 K = R / (R + Q)

 #filtered value
 x_{k|k} = x_{k|k-1} + K*(innovation)

 xprior <- update_xprior(x_{k|k})

 residuals <- update_residuals(measurement,x_{k|k})

 R = variance(residuals)
 Q = e^(-T)*R

 ...

I must be missing something...Also, I am not confident that what I am doing with T and the way I am calculating Q is correct...but I should be seeing smoothing regardless, since that is what they did in the document I am working off of.

So let's try to piece the signal model together first.

The state vector, $\mathbf{x}_k$ is just
$$
\mathbf{x}_k = [ x_{k-1}\ x_{k-2}\ x_{k-3} ]^T
$$
and the state update equation is just
$$
\mathbf{x}_k = \left[
\begin{array}{ccc}
3 & -3 & 1 \\
1 & 0 & 0 \\
0 & 1 & 0 \\
\end{array}
\right] \mathbf{x}_{k-1} +
\left[
\begin{array}{c}
1 \\
0 \\
0 \\
\end{array}
\right] v_k
$$
where $v_k$ is some driving noise which is Gaussian, zero mean, with variance $\sigma_v^2$.

Then the output equation is just
$$
x_k = \left[\begin{array}{ccc}
1 & 0 & 0
\end{array}
\right] \mathbf{x}_k + w_k
$$
where $w_k$ is Gaussian, zero mean, with variance $\sigma_w^2$.

If I set $\sigma_v$ to be $1$ and $\sigma_w$ to be 0.01 and generate one realization of the signal then I get the plot below (blue line)

This seems to be similar to the measurements you're graphing.

And the red line is what happens when I try to apply the Kalman filter to the signal model above.

Bottom line: I seem to get similar results to what you are seeing.

I suspect part of the problem is that the signal model is not stable, as the state transition matrix has eigenvalues:

  1.0000055+0.0000000i 0.9999973+0.0000048i 0.9999973-0.0000048i

the first of which is definitely outside the unit circle.

R Code Below

#27879

T <- 100

#H = [3 -3 1]
H <- matrix(c(1,0,0),c(1,3))

#F is just a tapped delay line
alpha <- matrix(c(0.9,0,0,0,0.9,0,0,0,0.9),3,3)
A <- matrix(c(3,1,0,-3,0,1,1,0,0),3,3)
B <- matrix(c(1,0,0),c(3,1))

sigma_v <- 1
v <- rnorm(T,0,sigma_v)
sigma_w <- 0.01
w <- rnorm(T,0,sigma_w)

x <- matrix(c(0,0,0),c(3,1))
z <- rep(0,T)

for (t in 1:T)
{
  z[t] <- H %*% x + w[t]
  x <- A %*% x +v[t]
}
par(mfrow=c(1,1))
plot(z,type="l", col="blue", lwd=5)

Q <- matrix(c(0,0,0,0,0,0,0,0,sigma_v^2),3,3)
R <- sigma_w^2

library("MASS") # For pseudo inverse ginv()

xkm1km1 <- matrix(rep(0,3*T+3),3,T+1)
xkkm1 <- matrix(rep(0,3*T),3,T)
K <- matrix(rep(0,3*T),3,T)
Pkm1km1 <- matrix(c(1000,0,0 ,0,1000,0, 0,0,1000),3,3)
zhat <- matrix(rep(0,T),c(T,1))

for (k in 1:T)
{
  xkkm1[,k] <- A %*% xkm1km1[,k]
  Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
  K[,k] <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  xkm1km1[,k+1] <- xkkm1[,k] + K[,k] %*% (z[k] - H %*% xkkm1[,k])
  Pkm1km1 <- (matrix(c(1,0,0,0,1,0,0,0,1),3,3) - K[,k] %*% H) %*% Pkkm1
  zhat[k] <- as.numeric(H %*% xkkm1[,k])
}
lines(zhat, col="red")


Output of IIR filter not matching the frequency response - Proximity to fs/4 After not having dealt with digital filters for a long time, I have now been playing around with filter design in octave and I am observing a behavior which I do not understand.
After designing a low-pass Chebychev filter with cheby1(1,3,0.4) and plot its frequency response, it can be seen that the -3dB point is exactly at the specified cutoff frequency.
However, when I test the filter with an input wave x=cos(2*pi*cutoff*t), using the filter function y=filter(b,a,x), the amplitude of the output signal can get significantly lower than the expected $0.707$, $0.63$ for this example in particular.
It seems that the output of the filter better approximates the plotted frequency response for cutoff frequencies either close to Nyquist frequency, or close to $0$, but as it gets closer to the middle point between these two, the output steps away from the expected attenuation.
What is the reason behind this?

From the octave documentation on cheby1.m you can see that the cut-off frequency in radians is given by pi*Wc, where Wc is the cut-off frequency input argument of cheby1.m. So for a function call [b,a] = cheby1(n,Rp,Wc), you should test the filter's gain at the cut-off frequency with an input sequence cos(Wc*pi*n). Furthermore, choose a sufficiently long input signal so that you can observe the filter's steady-state behavior.
Apart from that it's a matter of sampling the time domain function. Note that the samples will generally not occur at the maxima of the sinusoid. I.e., the maximum magnitude of the output signal
$$y[n]=|H(\omega_c)|\cos(n\omega_c+\phi(\omega_c))$$
will not necessarily attain its theoretical maximum value $|H(\omega_c)|$.


Separation of overlapping frequencies I have a signal with multiple frequencies, and two of them, one of which is my main frequency, overlap.

Are there any techniques that could separate two frequencies that almost overlap?
I can increase my sampling rate but not my window size

I have already tried SVD,PCA,ICA but none of the methods managed to decompose those two frequencies.

N=200
x = np.zeros(N)
y = np.zeros(N)
z = np.zeros(N)
k = np.zeros(N)
x[0] = np.cos(2*np.pi)*np.cos(2*np.pi*0.1)
y[0] = np.cos(2*np.pi)*np.sin(2*np.pi*0.1)
z[0] = np.sin(2*np.pi)*np.cos(2*np.pi*0.1)
k[0] = np.sin(2*np.pi)*np.sin(2*np.pi*0.1)
for i in range(1,N):
  x0=x[i-1]
  y0=y[i-1]
  z0=z[i-1]
  k0=k[i-1]
  x[i] = x0*np.cos(2*np.pi*0.31) - y0*np.sin(2*np.pi*0.31)
  y[i] = x0* np.sin(2*np.pi*0.31)+ y0*np.cos(2*np.pi*0.31)
  z[i] = z0* np.cos(2*np.pi*0.32)- k0*np.sin(2*np.pi*0.32)
  k[i] = z0* np.sin(2*np.pi*0.32)+ k0*np.cos(2*np.pi*0.32)
  y[i]+=0.3*(x[i]**2-z[i]**2)
  k[i]+=0.3*(x[i]*z[i])
jet= plt.get_cmap('jet')
colors = iter(jet(np.linspace(0,1,30)))
plt.figure(0)
w=30
for i in range(1,30):
  fourier = np.fft.rfft(x[i:i+w])
  freqs =np.fft.rfftfreq(len(x[i:i+w]))
  color=next(colors)
  plt.plot(freqs,abs(fourier), c=color)
plt.figure(1)
fourier = np.fft.rfft(x)
freqs =np.fft.rfftfreq(len(x))
plt.plot(freqs,abs(fourier), c='k')
plt.show()

In case your signal is fully stationary, just go ahead and use a "sinusoidal fitting" and take out the frequency that is not off your interest.
Then you have a disturbance free signal. Although since these signals are, I guess, not multiple integers of each other, your main frequency will also be low pass filtered. You can also use a notch filter, which will have similar behavior.

Please update us if you found a good answer.


Detect human speech in realtime audio on mobile phones I am looking to develop an Android app. As part of the functionality the app would require to randomly sample 3-5 seconds of audio and classify it as containing human speech or not. I understand that this concept is called Voice Activity Detection?

What would be the best way to implement this on a Mobile phone. I developed a basic system using energy based features and thresholds. I am hoping to find something less susceptible to noise, probably using features such as MFCC or formants? I did go through a number of papers, but most of them would require me to collect data and train models. Is there any library or framework I could use which would work in realtime?

I believe that speex at http://www.speex.org/ open source code has VAD inside. Try to see if you can see it and get some implementation ideas, with obaying their license.