Differentiation of sine in Fourier domain The derivative of $\sin(\omega_o t)$ is $\cos(\omega_o t)$.

The Fourier transform of $\sin(\omega_o t)$ is $\frac{\pi}{j}[\delta(\omega-\omega_o) - \delta(\omega+\omega_o)]$.

Differentiation in the time domain is equivalent to multiplying the transform by $j\omega$.

The transform of $\cos(\omega_o t)$ is $\pi[\delta(\omega-\omega_o) + \delta(\omega+\omega_o)]$.

What I don't understand is how multiplying the transform of $\sin(\omega_o t)$ by $j\omega$ gives you the transform of $\cos(\omega_o t)$.  I see how the $j$'s will cancel out, but how does the sign of that impulse get flipped?

I assume you mean the derivative with respect to $t$. In that case, the derivative of $\sin(\omega_0t)$ is not $\cos(\omega_0t)$ but $\omega_0\cos(\omega_0t)$. And luckily, this is also obtained via the Fourier transform relation you mentioned in your question:

$$\begin{align}\mathcal{F}\left\{\frac{d}{dt}\sin(\omega_0t)\right\}&=j\omega\cdot \mathcal{F}\left\{\sin(\omega_0t)\right\}\\&=\pi\omega[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]\\&=\pi[\omega_0\delta(\omega-\omega_0)-(-\omega_0)\delta(\omega+\omega_0)]\\&=\pi\omega_0[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]\\&=\omega_0\mathcal{F}\{\cos(\omega_0t)\}\end{align}$$

where I've used the fact that $f(\omega)\delta(\omega-\omega_0)=f(\omega_0)\delta(\omega-\omega_0)$ for any function $f(\omega)$ that is continuous at $\omega=\omega_0$. Consequently you have $\omega\delta(\omega+\omega_0)=-\omega_0\delta(\omega+\omega_0)$.


Auto-correlation of C/A code with 30dB noise I've a C# program that generates a GPS C/A code consisting of -1s and +1s.
Auto correlating the code with rotated versions of itself shows that the peak is 1023 and the next maximum value is 63 which agrees with what I expect.

Now my understanding is that in the GPS C/A signal the SNR is -30dB and the receiver is still able to detect a strong auto-correlation peak.

So I'm adding a random numbers between 0 and 30 (roughly 30dB) to one of the signals before the auto-correlation and now there's no strong peaks.

Where is the mistake in my reasoning?

using System;
using System.Collections.Generic;
using System.Linq;

namespace CrosssCorrelationTest
{
    class Program
    {
        static void Main(string[] args)
        {
            List<int> caCode = GenerateCaCode();
            var shiftedCaCode = new Queue<int>(caCode);
            AddNoise(caCode);

            double correlation = 0;
            double max1 = 0, max2 = 0;

            for(int i=1; i<1023; i++)
            {
                correlation = CrossCorrelate(caCode, shiftedCaCode.ToList());
                if(correlation > max1)
                {
                    max2 = max1;
                    max1 = correlation;
                }
                else if(correlation > max2)
                {
                    max2 = correlation;
                }
                RotateBy1(shiftedCaCode);
            }

            Console.WriteLine($"max1={max1}; max2={max2}");

            Console.ReadKey(true);
        }

        private static void AddNoise(List<int> caCode)
        {
            var r = new Random();
            for(int i=0; i<caCode.Count; i++)
            {
                caCode[i] += r.Next(30);
            }
        }

        private static double CrossCorrelate(List<int> caCode, List<int> caCode1)
        {
            int sum = 0;
            for(int i =0; i<caCode.Count; i++)
            {
                sum += caCode[i] * caCode1[i];
            }

            return sum;
        }

        private static void RotateBy1(Queue<int> shiftedCaCode)
        {
            int chip = shiftedCaCode.Dequeue();
            shiftedCaCode.Enqueue(chip);
        }

        private static List<int> GenerateCaCode()
        {
            int g1 = 1023; //10b'1111111111;
            int g2 = 1023; //10b'1111111111;

            var caCode = new List<int>();

            for (int i=1; i<=1023; i++)
            {
                int chip = (g1 ^ ((g2 >> 4) ^ (g2 >> 8))) & 1;

                caCode.Add(chip == 1 ? 1 : -1);

                g1 = (g1 >> 1) | (((g1 << 2) ^ (g1 << 9)) & 512);
                g2 = (g2 >> 1) | (((g2 << 9) ^ (g2 << 8) ^ (g2 << 7) ^ (g2 << 5) ^ (g2 << 2) ^ (g2 << 1)) & 512);
            }

            return caCode;
        }
    }
}

GPS has a processing gain of 43 dB assuming you correlate over 20 consecutive symbols. (GPS sends 20 C/A code sequences for 1 data bit: Chip rate 1.023 MHz, data rate 50 bps).  It looks like you are correlating over just one symbol which is at 1 KHz in which case the processing gain would only be $10Log_{10}(1023)= 30 \text{ dB}$.

However I don't think your noise is 30 dB stronger than your pre-correlation signal. The standard deviation for a uniform distribution with a range R is $\sigma = R/ \sqrt{12}$. If your range is 30, then the standard deviation would be 8.66. The DC offset does not effect standard deviation. The standard deviation for the C/A code that is +1/-1 is 1, so the predicted pre-correlation SNR would be 18.75 dB. In which case we should see peaks given an expected post-correlation SNR of 11.25 dB when correlating over 1 ms (1023 chips).

Here is the result I got from Matlab with the same test parameters (uniform noise 0 to 30, 1023 chips GPS C/A code +1/-1, correlated over one 1023 chip PRN sequence), where I also shifted 0 lag to be in the center of the plot.

Comparing my plot to yours I think your results are fine: I see a larger peak on the left side of your plot that is reasonably the expected correlation peak: Given the mean of the added noise is 15, the correlation and noise will be shifted up by 15, such that the mean of the peak would be 1023+15 = 1038. The peaks will vary by the post-correlation noise, which would have a standard deviation of  $8.66\sqrt{1023}= 277$ and be now Gaussian distributed due to the law of large numbers. Therefore 68% of the peaks would be expected to be in the range of 1038 +/- 277.

Cross Correlation of SV6 with SV6 + Noise over 1 ms

Added Noise

To Ben's point in the comments, this is a LOT easier to prototype from Matlab/Octave or Python:

Matlab/Octave Code:
(C/A Code generator cacode.m and further details on its implementation are available on this post GPS Coarse Acquisition PRN Codes)

s = cacode(6);  % C/A for SV6
s1 = s*2 -1;    % map to +1/-1
noise = rand(1,1023)*30;    % uniform noise
% circular cross correlation:
xcorr = ifft(fft(s1).*conj(fft(s1+noise)));
figure;
plot(real(xcorr));


signal sampling for 2nd time of sampled signal A sinusoidal signal of 600hz is sampled with 1khz, If the sampled signal is applied to an ideal low pass filter of 500 hz, what is the output?

thirI believe it is a nonsensical question, or at least a tricky, meant  to detect students who have  a shaky knowledge, and whose anwser deserves a lot of care.

First trick: if the $600$ Hz signal is sampled at  $1000$ Hz, it will endure aliasing (a proper sampling frequency would be above  $2 \times 600$ Hz). If the sampling process is perfect, it will appear as a sine at the mirror frequency $2*\frac{1000}{2}-600 = 400$ Hz.

The  signal (that is assume real)) is now  in the $[0,500] $ Hz band. An ideal-low pass cutting at $500$ Hz would thus be an all-pass, second trick, and would be utterly  hard to design, numerically, in that case (third trick).

So in theory, the  output signal would be the same (aliased) signal as the one from the first stage, fourth trick in a way.


Difference between speech and music signal What are the difference between speech and music signal? Can anyone help?

I am doing project in audio processing. In that project I want make decision, that is which type of signal is given as input whether it is speech signal or music signal. I searched a lot in internet. But i could not find the required details.

If anyone knows about this please share the difference between the speech and music signal.

Not really an answer, but usually when people speak there are pauses between words, sentences, thoughts, etc. Music (with multiple instruments) is fairly continuous. Even if one instruments stops another is probably still playing. So there would be less lulls in the audio stream. Just a guess though.

Another thing you could do is look at frequencies, voices are typically in the 300Hz-4kHz range (this is what telephones originally used as bandwidth) whereas musical instruments are usually much broader. A frequency outside the human range could indicate instruments  I randomly found that image here


Matlab - How to create an autocorrelogram using a spike-train I have already posted this question on MATLAB's question forum.  However, I haven't received any responses, and am hoping to find some folks who may be a bit more signal processing savvy on here.  Here goes... I'm wondering if someone could help clear up some confusion I am having surrounding the xcorr() function since I am struggling to parse how it works exactly.

Here is the matlab description:

  r = xcorr(x) returns the autocorrelation sequence of x.

I am trying to make an autocorrelogram out of spike train data. Thus, I am inputting timestamps into xcorr() like so... xcorr(spiketimes). I am operating under the assumption that autocorrelating time-stamps for events (spikes) should lead to an autocorrelogram representative of the the autocorrelation between interspike intervals. This autocorrelogram, I believe, should provide the same info as an autocorrelogram of 0's and 1's (0 = no spike, 1 = spike) inputted into a vector who's length is the size of the recording (in milliseconds).  When I run that line of code, however, I find that there's a peak towards the length of the timestamps vector. This doesn't make sense. The peak should be at lag zero (where the vector is 100% correlated with itself). What am I doing wrong here? Or is my understanding of an autocorrelogram incorrect? Does anyone here have any experience using xcorr(), and have any advice?

Also, this is less important and a bit technical, but for fellow matlab users, I'm not totally sure how to adjust bin width. If anyone has some pointers, I'd appreciate it. Finally, the fact that I can't center the plot that's created from xcorr() at zero is driving me crazy. How can I fix this?

When I run that line of code, however, I find that there's a peak towards the length of the timestamps vector.

So I'm assuming you're running:

r = xcorr(x);

and you're perplexed by the peak at r(length(x)) ?  That's because xcorr returns a vector that is the autocorrelation from a lag of -length(x) to a lag of +length(x) (unless told otherwise) so that the middle value is the zero lag.

However you seem interested in finding correlations between the pulse times.  That's a whole other kettle of fish... especially when there are pulse trains of different periods or there are missing data. To see the complexity that can happen, check out this paper.

What are you trying to achieve?

To see that there's a big difference between the autocorrelation function of the pulse arrival times and the actual pulses themselves, I've generated an example in R.  The resulting autocorrelation functions are plotted below.

--

Ok, so you have a periodic(ish) pulse train and you want to figure out what the pulse period is.

The autocorrelation of the pulse times will not help. This is because the times are just a sorted array, so you'll just get a straight line as the autocorrelation.
The autocorrelation of the pulses themselves should help.  However, you'll get a peak at 0 lag, and at lags of $nP$ where $P$ is your pulse period and $n$ is an integer.
The post you reference is looking at the histogram of the inter-pulse times (time differences between successive pulses).  Depending on how varied your pulse train is, this might be the best approach.

If I change my simulation to use a pulse that is periodic (with zero noise on the pulse times!), then I get the second lot of code below, and the results are graphed here:

For this particular run, the Tperiod = 9.350267.  As you can see, both the ACF of the pulses and the histogram method get a peak close to where we want it (though quantization effects mean we don't get exactly at 9.350267; better selection of "bins" might improve that).

R Code Below

#26547

Npulses <- 100
Tmax <- 1000

pulseTimes <- sort(runif(Npulses,0,Tmax))

# install.packages("Hmisc")
library(Hmisc)

xc_times = acf(pulseTimes, plot = FALSE)

times <- 0:.1:Tmax

pulses <- times*0

for (pulseNo in 1:length(pulseTimes))
{
  idx <- which.min(abs(pulseTimes[pulseNo] - times))
  pulses[idx] = 1
}

xc_pulses = acf(pulses, plot = FALSE)

plot(xc_times$lag,xc_times$acf, type="l", ylim=c(-0.1,1))
lines(xc_pulses$lag, xc_pulses$acf, col="red")
legend(15, 0.9, c("ACF of times", "ACF of pulses"), lwd=c(2.5,2.5),col=c("black","red"))

Second version of R code below

#26547

Npulses <- 100
Tmax <- 1000

Tperiod <- runif(1,9,10)
Tstart <- runif(1,0,1)

#pulseTimes <- sort(runif(Npulses,0,Tmax))
pulseTimes <- Tstart + (0:Npulses)*Tperiod

# install.packages("Hmisc")
library(Hmisc)

xc_times = acf(pulseTimes, plot = FALSE)

times <- 0:.1:Tmax

pulses <- times*0

for (pulseNo in 1:length(pulseTimes))
{
  idx <- which.min(abs(pulseTimes[pulseNo] - times))
  pulses[idx] = 1
}

xc_pulses <- acf(pulses, plot = FALSE)

histogram <-  hist(diff(pulseTimes), breaks=xc_pulses$lag, plot=FALSE)

plot(xc_times$lag,xc_times$acf, type="l", ylim=c(-0.1,1))
lines(xc_pulses$lag, xc_pulses$acf, col="red")
lines(histogram$mids, histogram$counts / max(histogram$counts), col="green")

legend(12, 1.0, c("ACF of times", "ACF of pulses", "Histogram of inter-pulse times"), lwd=c(2.5,2.5, 2.5),col=c("black","red", "green"))


adding phase noise in signal thorugh freuquency I want to implement a simple signal for phase noise. According to my supervisor, i shoud be adding phase noise in the frquency. This is my signal. A = cos(2*pi*fc*t) and signal with phase noise B = cos(2*pi*(fc+deltaf)*t). here how shoud i define delta f?

After that when i analyse both signals using fft, spectrum of B signal should be broader than spectrum of A signal.

First - this will depend on the nature of the phase noise you wish to model. If deltaf is a constant then it is just a frequency shift and really isn't considered phase noise. What you really need is deltaf to be a noise vector the same length as t. In matlab format you can do this by:

t=0:255;
fc=0.2; %just an example
deltaf = 0.001*randn(size(t));
B=cos(2*pi*(fc+deltaf).*t);

Notice that the frequency will have a finite variance in this case - but the phase variance will increase linearly over time - this may or may not be what you want. To have a finite phase variance you would do something like

t=0:255;
fc=0.2; %just an example
ph_n = 0.0001*randn(size(t));
B=cos(2*pi*fc*t+2*pi*ph_n);

Typically you can analyze these signals by using trigonometric identities and then assuming the noise portion is small or by using Taylor series expansions for the noise components.


Sampling multiple sensors at different rates I have a question regarding how best to sample data from multiple sensors when they are all updating at different/changing rates.  I am new to data acquisition so please bear with me.  For my Master's research, I am trying to build a dynamometer to measure the performance of a belt transmission used with the Baja SAE team at my university.  This is a VERY low budget project, so I am trying to collect data with arduino/raspberry pi if possible.  The test setup consists of a drive pulley attached to an engine and a driven pulley with a friction brake behind it.  I am trying to measure the RPMs for both shafts, and the torque on the driven shaft.  To measure RPM's, I have magnets mounted on the shaft and I am using hall effect sensors to measure the time between rotations.  To measure the torque I have strain gages with slip rings.

The problem is that all sensors update at different rates.  The drive pulley rotates at about 60Hz, and the driven can vary from 20-120Hz, with the values being measured once per rotation.  The A/D converter for the strain gages reads at a constant 125Hz (the input voltage has a low pass filter to remove engine and other higher frequency noise-I am only interested in time averages).  My goal is to be able to get a plot of the steady-state driven pulley torque vs speed.

I tried initially reading all of the sensors simultaneously with a raspberry pi, but all of the sensors are read with interrupts in the python code and the timing was off since it is not a real time device (the RPM data was VERY noisy).  If I read each RPM sensor and the A/D converter individually with an arduino, it has no problem keeping up and the data looks clean.  Now that I have all of the sensors reading, I am trying to put it all together.  Since the reaction time of the transmission to a changing load is fairly slow, I think I could capture enough information if I sampled at 20-25Hz.

My question is if it would work to have separate arduinos reading each of the sensors continually and updating global variables in their RAM, then calling these current running values from the raspberry pi at a fixed sampling rate (20-25Hz).  I didn't know if it would be a problem that I don't have timestamps for when the global variables last updated.  I will not be using any kind of frequency analysis on the data.  Also, to reduce noise from the strain gage A/D converter, could I have that global variable be a running average of the last, say 5, data points so that I am capturing information from each sample even though I am only calling values at 1/5 of its sample rate?

My suggestion is to use a single Arduino for all the sensors, create timestamps for the data and transmit it via I2C to whatever sink you want. Your sampling rates are so low, that a single uC should do the job.
Then do a generous oversampling to make data from different sensors easy to synchronize.


Calculating an output of a system (Z- transform question) I have a following question to answer:

An LTI system is described by its impulse response h[n].  For input x[n] it gives output y[n].

$$h[n] = u(n) - u(n-N) $$
$$x[n] = u(n) - u(n-M)$$

I want to calculate it's output function. I have used the Z-transform property of convolution in time being multiplication in Z-domain to calculate Y(z) which looks like this $$Y(z) = \frac{1-z^{-M}}{1-z^{-N}} $$
However I have trouble to apply inverse Z-transform to it. Is there some nice way I can transform this resulting function to time(sample) domain?

For this specific example, using a partial fraction expansion (PFE) would be an overkill, so yo should better consider using the formula :

$$ 1 - z^{-N} = (1- z^{-1}) (1 + z^{-1} + z^{-2} + ... + z^{-N+1} $$

to simplify the Z-transform $Y(z)$ of your output as:

$$\begin{align}
Y(z) &= H(z) X(z) = \frac{1 - z^{-N} } {1-z^{-1} } \frac{ 1 - z^{-M} }{ 1-z^{-1} } \\ \\
&= \left( 1 + z^{-1} + z^{-2} + ... + z^{1-N} \right) \left( 1 + z^{-1} + z^{-2} + ... + z^{1-M} \right)\\
Y(z) &= 1 + c_1 z^{-1} + c_2 z^{-2} + ... + c_{N+M-2} z^{2-N-M} \\
Y(z) &= \sum_{n=0}^{K} y[n] z^{-n} = 1 + c_1 z^{-1} + c_2 z^{-2} + ... + c_{N+M-2} z^{2-N-M} \\
\end{align}$$

Now as you can see, $Y(z)$ is a simple polynomial in $z^{-1}$. You can either explicitly multiply those two polynomials in parenthesis to find those coefficients $c_k$ to get th sample values $y[n]$ as indicated by: $$y[k] = c_k ~~~,~~~ k = 0,1,2,...,N+M-2$$

you can easily get those (coefficient $c_k$) values by the following convolution:
$$ y[n] = c[n] = a[n] \star b[n]$$ where $a[n]$ and $b[n]$ are all ones of length $N$ and $M$ respectively.

Now if you insist on a PFE method, the you could also do that as follows:
$$\begin{align}
Y(z) &= H(z) X(z) = \frac{1 - z^{-N} } {1-z^{-1} } \frac{ 1 - z^{-M} }{ 1-z^{-1} } \\ \\
&= \frac{1 - z^{-N} - z^{-M} + z^{-N-M} } { (1-z^{-1})^2 }\\ \\
&= \frac{1} { (1-z^{-1})^2 } - \frac{z^{-N}} { (1-z^{-1})^2 } - \frac{ - z^{-M}  } { (1-z^{-1})^2 } + \frac{z^{-N-M} } { (1-z^{-1})^2 }\\ \\
\end{align}$$

Now as can be seen, I have separated the expression into four pieces each of which is a delayed version of the first one (on the left without a delay). So I will express $y[n]$ as a sum of those four delayed functions, denoting the first one as $f[n]$ we have :

$$y[n] = f[n] - f[n-N] - f[n-M] + f[n-M-N]$$ where the function $f[n]$ is the inverse Z- transform of
$$ f[n] = \mathcal{Z}^{-1} \{ \frac{1}{(1-z^{-1})^2} \} $$

looking from a table or simply finding it yourself, one can find that
$$f[n] = u[n] \star u[n]  = r[n]$$

where $r[n]$ is the following function :
$$r[n] = (n+1) u[n]$$

hence the result is simplified as:

$$ y[n] = r[n] - r[n-N] - r[n-M] + r[n-M-N] $$

instead of denoting the result in terms of ramp functions, you can also simplfity the shifted sum formula through its convolution operatos as:
$$ y[n] = (u[n] \star u[n]) -(u[n] \star u[n-N])-(u[n] \star u[n-M])+(u[n] \star u[n-N-M])$$
$$ y[n] = u[n] \star \left(  u[n] - u[n-N] - u[n-M] + u[n-N-M] \right)$$

denoting the finite length sequence on the right by $g[n]$, then the output can be found to be

$$y[n] = u[n] \star g[n] = \sum_{k=-\infty}^{n} g[k] = \sum_{k=0}^{n} g[k]  $$

The rightmost sum occurs assuming $g[n]$ is a causal sequence.


Auto-correlation of the sum of two generic signals Be $x[n]$ and $y[n]$ two generic discrete-time signals. Given $s[n] = x[n] + y[n]$ I want to evaluate its autocorrelation $R_s[l]$.

By definition (https://en.wikipedia.org/wiki/Cross-correlation):
$$R_s[l] = \sum^{+\infty}_{m = -\infty}s^*[m]s[m+l]$$
Which expands to:
$$\dots = \sum^{+\infty}_{m = -\infty}(x[m]+y[m])^*\,(x[m+l]+y[m+l])$$
Since $(a + b)^* = a^* + b^*$ :
$$\dots = \sum^{+\infty}_{m = -\infty}x^*[m]x[m+l]+y^*[m]y[m+l]+x^*
[m]y[m+l]+y^*[m]x[m+l]$$
Using the definition again:
$$\dots = R_x[l] + R_y[l] + R_{xy}[l] + R_{yx}[l]$$

Is this result correct?

This is correct. You can use $R_{xy}[l]=R_{yx}^{*} [l]$ to simplify a little further to

$$R_{ss}[l] = R_{xx}[l]+R_{yy}[l]+2 \cdot  \Re\{ R_{xy}[l] \} $$


Why when the SINR of transmitting signal is very large ,the device cannot decode the received signal at the same time? $A$ transmit $x$ and receive the $y$ at the same time.If the SINR for transmitting $x$ is very large,$A$ cannot decode the $y$ signal.Why?

In my opinion,If the SINR for transmitting $x$ is very large,it means that SINR for receiving $y$ is very small,so the error probability for receiving $y$ will be very large,so the $A$ cannot decode the $y$ signal.

Is my thinking right?

A transmit x and receive the y at the same time,so of course it is a full duplex,i think your thinking is right.

If two antennas,A & B, are very close,and A is for transmitting,B is for receiving.B can receive a very strong signal from A,that is,if B is not going to to receive the signal from A,but C.The SINR for receiving C will be very small. I mean

$SNR_{AinB}=\frac{A \ power}{C \ power +noise}$,if A and B are very close,the $A \ power$ will become very large,that is,$A \ power >> C \ power$

And if B want to receive the signal from C in fact,the SINR will become

$SNR_{CinB}=\frac{C \ power}{A \ power +noise}$,and as i just mentioned, $A \ power >> C \ power$,so we can know $SNR_{CinB}$ is very small.

SNR is defined as the ratio of signal power to the noise power,so i think it can also apply on transmitting,not only on receiving.Nevertheless,most people will apply on receiving

If SNR is very small,the information cannot transmit to the receiver correctly,so of course we cannot know,or decode, the real information


Discrete Time Fourier Analysis Suppose we're given the following:
$ x[n] = 2 + (-1)^n $, and are given the impulse response $ h[n] = u[n] a^n $, of an LTI system where $ |a| < 1$. We're asked to find the output $y[n]$, if $x[n]$ is input into $y[n]$.

My approach to this problem was to find the frequency response to $x[n]$ and the frequency response to $h[n]$, multiply the two, and then find the inverse fourier transform. However, this method was quite complicated, and I ended up not being able to solve it, since I essentially got back to a form where it was just $x[n]$ convolved with $h[n]$, which is difficult to compute by hand.

The answer key has a far simpler approach where $x[n]$ is written as $2 e^{j 0 n} + e^{ j \pi n}$, evaluates the frequency response $$ H = \frac{1}{1-a e^{-j \omega}} $$ at $\omega$; values of $0$ and $\pi$, multiplies and adds the result to yield:

$$ y[n] = \frac{2}{1-a} + \frac{1}{1-a} (-1)^n $$

Why can this be done? I'm not sure I follow why $y[n]$ can essentially by written as $H(\omega)$ evaluated at a value of $\omega  * x[n]$.

Can this always be done for LTI systems--even continuous systems?

The sequences $e^{j\omega_0n}$ are eigensequences of discrete-time LTI systems, i.e., the response to such a sequence is the same sequence scaled by a complex constant (the eigenvalue). This can be shown as follows:

$$\begin{align}y[n]&=(h\star x)[n]\\&=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\\&=\sum_{k=-\infty}^{\infty}h[k]e^{j\omega_0(n-k)}\\&=e^{j\omega_0n}\sum_{k=-\infty}^{\infty}h[k]e^{-j\omega_0k}\\&=e^{j\omega_0n}H(\omega_0)\end{align}$$

where $h[n]$ is the system's impulse response, and $H(\omega)$ is its frequency response. A completely analogous proof can be given for continuous-time systems.


Could there be any reason to prefer convolution-based calculation of autocorrelation? Theoretically both of ways of calculating autocorrelation function are identical: strightforward convolution and Fourier-based method where we use FFT/iFFT in practice. And as it is well known, the computational complexity is $\mathcal O(N^2)$ and $\mathcal O(N\log N)$ for them respectively.

So the question is: are there any reasons why one could prefer convolution over the FFT-based method?

I could imagine only the memory-related argument (e.g. when we don't have any extra memory for calculations), but thinking of a regular PCs, at the time when memory constrains start to play role, the time complexity of $N^2$ will already be a killer for the first method.

Are there any particular numerical issue of the FFT that one needs to be aware of here (and choose convolution to avoid it)?

And as it is well known, the computational complexity is $\mathcal{O}(N^2)$ and $\mathcal{O}(N \log N)$ for them respectively.

$\mathcal{O}$ notation ignores any constants that determine exactly how fast those functions run. Depending on the constants out front, it may be faster to compute the convolution (or correlation) directly.

For example, the Fourier transform method to convolve may take $10^5 N \log N$ seconds while the direct method may take $10^{-3} N^2$ seconds. If $N = 100$, the direct method is faster.

In general for large $N$, the Fourier method tends to be faster with speedups of up. For more detail on computing which method is faster, see scipy PR #5608.

Below is a graph that times the Fourier transform method of convolution as fftconvolve and the direct method as convolve. We see that the direct method is faster for many arrays of practical size.

Note: this graph is when convolving 1D signals. For more higher dimensional signals (2D, 3D, etc) fftconvolve tends to be faster more often

To time these graphs, I used scipy.signal.fftconvolve and np.convolve with the default parameters. As of scipy version 0.19, scipy.signal.convolve will choose between the direct and Fourier transform methods to choose the faster method.


Understanding Asymptotic Equipartition Property I have some problems in understanding the precise meaning of the Asymptotic Equipartition Property, related to a large number n of independent and identically distributed random variables (X1, X2, ..., Xn). Here an example is shown (here you may find the complete material):

I do not understand the precise difference between the sentences:

"Clearly, it is not true that all 2n sequences of length n have the same
probability."
"One might be able, however, to predict the probability of the sequence that is actually observed. The question is: what is the probability p(X1,X2,...,Xn) of the outcomes X1,X2,..., Xn, where X1,X2,..., Xn are i.i.d. according to p(x)"
"the typical set has probability close to 1, all the elements of the
typical set have same probability"

Then I read this sentence, which increases the confusion in my mind.

The sentences you've mentioned in your question all have a quite different meaning. Let me try to explain them one by one:

  Clearly, it is not true that all $2^n$ sequences of length $n$ have the same probability.

This refers to the example of a sequence of binary i.i.d. random variables. If $p(1)\neq p(0)$ then it is clear that sequences with different numbers of ones and zeros must have different probabilities, i.e., not all sequences have the same probability.

  One might be able, however, to predict the probability of the sequence that is actually observed. The question is: what is the probability $p(X_1,X_2,...,X_n)$ of the outcomes $X_1,X_2,\ldots, X_n$, where $X_1,X_2,\ldots, X_n$ are i.i.d. according to $p(x)$?

This question is answered by the asymptotic equipartition property (AEP). It turns out that the typical sequences do have approximately the same probability as $n$ gets large, namely $2^{-nH}$.

  The typical set has probability close to 1, all the elements of the typical set have [the] same probability.

That's also what the AEP says. As mentioned before, the typical sequences have approximately the same probability, and the aggregate probability of the typical set approaches $1$ as $n$ gets large.

Finally, the statement

  The concept of a typical set is generally different from that of a high probability set.

is indeed confusing. What is meant is the fact that there can exist individual sequences that are not part of the typical set but which have a higher probability than the typical sequences. Still, the aggregate probability of the typical set approaches $1$ for increasing $n$, and there are too few atypical sequences with high probability to have any significance.


Question about Prony's method I want to use Prony's method to fit a signal given by the sum of $p$ complex exponentials

$x[n]=\displaystyle{\sum_{k=1}^p} A_k{\rm e}^{\,{\rm j}\theta_k} {\rm e}^{2\pi\, {\rm j} f_k (n-1)}$.

The problem is that Prony's standard formulation (Hildebrand, Introduction to Numerical Analysis, pp. 378-382) assumes that all the samples $x[n]$ with $n=1,\ldots,N$ are available, however in some applications $x[1]$ can't be measured. It seems that in this case Prony's method can't be applied. Am I missing something?

Drop the first sample - apply the formulation as if x[2] was x[1].


Identification of a given barcode as either 1 D or 2 D I am working on an algorithm (in MATLAB) to identify/classify a given bar code (image) as 1 D or 2 D. I have considered using the contrast of THE IMAGE matrix as a parameter or using the FFT of the image.

I would be grateful for any ideas or suggestions on what you consider the best approach for this problem

I think this is a fairly simple problem. Just compute the absolute gradient histogram (you, probably, don't even need to normalize). If it is concentrated around one angle, this should indicate a binary code, whereas more or less equally distributed histograms should be the opposite. Of course, this approach assumes that you always see barcode images, i.e. no negative class.

Ok, so I did a sample implementation for you, considering the above algorithm I described. Of course you could always make it better, but for the real-life images I tested, it gives reasonable results:

function [type] = classify_barcode_type(I, mingrad)

% convert to grayscale
if (size(I,3)>1)
    I = rgb2gray(I);
end

% smooth the image
I = imgaussfilt(I, 2);

% compute sobel gradients
[Gmag,Gdir] = imgradient(I);

% upside down are the same (180 deg. range)
gdir = mod(Gdir+180,180);

% threshold the gradients: reduces noise
gradSelect = find(Gmag>mingrad);
numGrad = length(gradSelect);

% histogram of gradients
[N, ~] = histcounts(gdir(gradSelect), 180);

% 'relMag' percent of the gradients are pointing towards a single dir.
relMag = max(N)./numGrad;

if (relMag<0.05)
    type = 'qr';
else
    type = 'linear';
end

end

Sample images for an Aztec code and Code 128 barcode are as follows:

You could use the MATLAB function as follows:

% read the image
img = imread('real6s.png');

% apply custom rotation for robustness test
imgr=imrotate(img,125);

% classify
[type] = classify_1d_barcode(imgr, 40)

Note that, even though rotation creates many artificial edges, the approach remains to be robust, even in presence of severe rotations. In real life, you won't have the synthetic edges appearing due to rotations.

The figures below are plots of normalized gradient histograms ($\frac{\mathbf{N}}{\lVert \mathbf{N} \rVert}$ - also in the code). Note the difference between the aztec (QR) code and the linear barcode:

Note that linear barcode is rotated by 90 deg and this is clearly visible in the normalized angle histogram plot. You can get the rotation angle this way. Even though this example was for 90 deg rotation, you will get the same / similar results with any other rotation - theoretically the same results actually. Threshold (or the decision boundary) is indicated by the green line.


What is "filter periodization"? A library defines periodize_filter_fourier, which is an equi-spaced averaging formulated by
$$
v_f[k] = \sum_{i=0}^{\text{n_periods}-1} h_f[i\cdot N + k],
$$
where $v_f$ is periodization of $h_f$, $N=\text{len}(h_f)/\text{n_periods}=\text{len}(v_f)$, $k=[0,1,...,N-1]$, and $1/\text{n_periods}$ seems missing.
It's applied on the frequency-domain (DFT) Morlet wavelet. What's the purpose, and what's "periodizing" to do with it?

Example with h = [1, 2, 3, ..., 12]:
# n_periods == 2
v_f = [[1, 2, 3, 4, 5, 6],
       [7, 8, 9, 10, 11, 12]].mean(axis=0)
    = [4, 5, 6, 7, 8, 9]

# n_periods == 3
v_f = [[1, 2, 3, 4],
       [5, 6, 7, 8],
       [9, 10, 11, 12]].mean(axis=0)
    = [5, 6, 7, 8]

TL;DR The idea is to modify the continuous-time representation of the filter, $f(t)$, such that it's suitable for discrete-time sampling - mainly by accounting for aliasing via "folding" in Fourier domain. The operation is also equivalent to (and used for) subsampling.
"Periodization" refers to the fact that sampling $f(t)$ at uniform intervals $s$ corresponds to making its Fourier transform $2\pi/s$-periodic.

Periodizing
A uniform sampling of $f$ with interval $s$ can be expressed as a weighted Dirac sum:
$$
f_d(t) = \sum_{n=-\infty}^{\infty} f(ns) \delta (t - ns) \tag{1}
$$
Fourier transform of $\delta(t - ns)$ is $e^{-jns\omega }$, so the Fourier transform of $f_d$ is
$$
\hat f_d (\omega) = \sum_{n=-\infty}^{\infty} f(ns) e^{-ins\omega} \tag{2}
$$
($f(ns)$ is the "weights" and thus constant in transform over $t$). It can be shown${}^1$ that this is same as
$$
\hat f_d (\omega) = \frac{1}{s}\sum_{k=-\infty}^{\infty} \hat f \left( \omega - \frac{2\pi}{s}k \right) \tag{3}
$$
That is, sampling $f$ at intervals $s$ is equivalent to making its Fourier transform $2\pi/s$-periodic by summing all its translations $\hat f (\omega - (2\pi / s)k)$. This makes more sense with what follows.

Sampling
(Shannon-Whittaker sampling theorem) If the support of $\hat f$ is bounded in $[-\pi/s, \pi/s]$, then
$$
f(t) = \sum_{n=-\infty}^{\infty} f(ns) \phi_s(t - ns), \tag {4}
$$
where
$$
\phi_s(t) = \frac{\sin{( (\pi/s)t )}}{(\pi / s) t}. \tag{5}
$$
Note, if $n \neq 0$, the support of $\hat f(\omega - (\pi / s)n)$ does not intersect the support of $\hat f (\omega)$ because $\hat f(\omega) = 0 $ for $|\omega| > \pi / s$. So $(3)$ implies
$$
\hat f_d (\omega) = \frac{\hat f(\omega)}{s}\ \ \text{if}\ \ |\omega| \leq \frac{\pi}{s} \tag{6}
$$
Fourier transform of $\phi_s$ is $\hat \phi_s = s \bf{1}_{[-\pi/s, \pi/s]}$ (indicator function, or "box"). Since the support of $\hat f$ is in $[-\pi/s, \pi/s]$, it follows from $(6)$ that $\hat f(\omega) = \hat \phi_s(\omega) \hat f_d(\omega)$. Its inverse Fourier transform gives
$$
\begin{align}
f(t) = \phi_s \star f_d(t)
&= \phi \star \sum_{n=-\infty}^{\infty} f(ns) \delta(t - ns) \\
&= \sum_{n=-\infty}^{\infty} f(ns) \phi_s(t - ns) \tag{7}
\end{align}
$$
This happens to be a proof of the theorem.

Sampling $\Leftrightarrow$ Periodizing
This can be seen as follows; suppose $\hat f$'s support is in $[-\pi/s, \pi/s]$. Then:

$(a)$ is simply $f(t)$ and its F.T. Sampling it at uniform interval $s$ has the effect of adding translations of its F.T. (i.e. $(3)$), which is what's shown - original $\hat f(\omega)$ translated and centered at multiples of $(2\pi /s)$. These "duplicates" extend out to infinity; we've thus "periodized" the Fourier transform of $f(t)$.
Per $(7)$, we recover $f(t)$ from $f_d(t)$ by convolving with $\phi_s$ (or multiplying in Freq domain):

Subsampling
Let's pose it as a problem: what effect does subsampling (i.e. taking every $k$-th sample) of a finite sequence $x[n] = [4, 7, ..., -1]$ have on its Fourier transform? Note this is equivalent to scaling the sampling interval: $s \rightarrow ks$. Suppose we wish to subsample by 2.
To answer this, we first visit aliasing.
Aliasing
In case $\hat f$ isn't bound to $[-\pi/s, \pi/s]$, the resulting $f_d(t)$ will alias. From $(3)$ and proof of $(4)$, we can predict this manifests as overlaps in Fourier domain:

The translations are still spaced by $(2\pi / s)$, while support of $\hat f$ exceeds it, yielding overlaps. The sampled $f_d(t)$ will then no longer resemble $f(t)$:

Subsampling (cont'd)
If $k=2$, we double our $s$, and thus halve the support $[-\pi/s, \pi/s]$. To subsample without aliasing, $\hat f$'s support must thus lie in $[-\pi/(2s), \pi/(2s)]$ (per original $s$):

Now suppose it isn't, and the support lies exactly in $[-\pi/s, \pi/s]$. Then:

What is the resulting $\hat f_d$? Well, it's simply the negative half of $\hat f(\omega)$ added to its positive half, and vice versa ("folding") - then divided by 2 (s *= 2 in eq $(3)$).
This is exactly what
$$
v_f[k] = \frac{1}{k} \sum_{i=0}^{\text{n_periods}-1} h_f[i\cdot N + k], \tag{8}
$$
(expression in question) is doing, but instead with DFT bins, where the de-facto support w.r.t. continuous-time input is $[-\pi/s, \pi/s]$ (with or without aliasing).
Put one way, subsampling simulates the effects of aliasing. In code, the following  holds (assume fold() implements $(8)$):
ifft(fold(fft(x), k)) == x[::k]

Periodizing a filter
Does $(8)$ serve any purpose beyond subsampling? Correctness: if we disregard folding and pretend aliasing isn't there, result's even worse than with aliasing.
More generally, given some arbitrary continuous-time function of a filter in Fourier domain, $\hat \psi(\omega)$, how do we sample it such that the entire support (assuming finite) is included, and that we use exactly $N$ samples? One approach is, "sample lots, then subsample" - which is exactly what Kymatio's morlet_1d does. (The alternative is to determine correct sampling bounds and get it right the first time.)
Once we have the correct $\hat \psi_d(\omega)$ that accounts for aliasing and captures at least one full period of the continuous-time function (i.e. $[-\pi/s, \pi/s]$), then $\hat \psi_d(\omega)$ is the correct continuous-time function, whether without or with aliasing, and taking IFFT we get closest to what we'd get with sampling the continuous-time $\psi (t)$ with rate $f_s=1/s$ (over one continuous-time period).
References
Figures and formulae (with much accompanying wording copied verbatim) taken from S. Mallat's Wavelet Tour, sections 3.1.0-3.1.2, where further explanations and proofs are given (including to footnote 1).


Correlation of a signal I have one sample for a signal. This sample is a vector of length 384. I need to calculate the correlation matrix for this signal,So I need many samples for the same signal.
How can i generate these samples form the given sample using Matlab?

The Toeplitz matrix is used to compute correlation and convolution using matrix multiplication. Below is a graphic showing how to use a Toeplitz matrix specifically to perform convolution using matrix multiplication. A cross-correlation function can also be done following the same process by simply reversing the order of the coefficients in the multiplied vector (in this case we would multiply by [h3 h2 h1] to perform correlation).

This is clearer by observing the similarity and differences in the formulas for discrete time convolution and cross-correlation:

CONVOLUTION:

$$r[n]*h[h] = \sum_{m=-\infty}^\infty r[m]h[n-m]$$

CROSS CORRELATION:

$$\rho_{rh}[n] = \sum_{m=-\infty}^\infty r[m]h[m+n]$$

The code in the blue cell in the graphic is the implementation in either Matlab or Octave. This applies to auto-correlation or cross-correlation as an alternate to using the xcorr() function directly.

The utility of doing this is clearer in the following post, where I used such an operation to help demonstrate the underlying operation of the LMS equalizer:

Compensating Loudspeaker frequency response in an audio signal


can $e^{j\phi_0}$ be incorporated into the sine and cosine terms? This is from "communication system" by Carlson fifth edition page 109:
If transfer function of a channel be like this:
$H(f)= Ae^{j\phi_0}e^{-j2\pi ft_g} $
and input to this system be:
$x(t)=x_1(t)\cos(\omega_c t) - x_2(t)\sin(\omega_c t)  $
will output be like this?
$y(t) = Ax_1(t-t_g)\cos[\omega_c (t-t_g)+\phi_0] - Ax_2(t-t_g)\sin[\omega_c (t-t_g)+\phi_0]$
my problem is with this sentence in the book:

$e^{j\phi_0}$ can be incorporated into the sine and cosine terms...

Edit:
picture from book:

Your problem might be the unmentioned fact that the given frequency response $(4)$ in Carlson's book is only valid for positive frequencies. The complete definition should be
$$H(f)=\begin{cases}Ae^{j\phi_0}e^{-j2\pi ft_g},&f>0\\Ae^{-j\phi_0}e^{-j2\pi ft_g},&f<0\end{cases}\tag{1}$$
Now we have $H(f)=H^*(-f)$, which is necessary for the corresponding system (channel) to be real-valued.
If you split the system into two subsystems $H_1(f)=Ae^{j\textrm{sgn}(f)\phi_0}$ and $H_2(f)=e^{-j2\pi ft_g}$, where $H_1(f)$ is a (scaled) phase shifter, and $H_2(f)$ is a pure delay, it's easy to see the impact of the total system on the input signal.
Consider the input signal $x_1(t)\cos(\omega_ct)$. The output of the scaled phase shifter is $Ax_1(t)\cos(\omega_ct+\phi_0)$, if the signal $x_1(t)$ is a lowpass signal with a maximum frequency smaller than the carrier frequency $\omega_c$. Delaying that signal by $t_g$ gives
$$y_1(t)=Ax_1(t-t_g)\cos[\omega_c(t-t_g)+\phi_0]\tag{2}$$
The result for the signal $x_2(t)\sin(\omega_ct)$ follows in a completely analogous way.
In sum, the two necessary conditions for the given formula to be true are

$H(f)$ corresponds to a real-valued system, i.e., $H(f)=H^*(-f)$, and
$x_1(t)$ and $x_2(t)$ are lowpass signals with maximum frequencies smaller than the carrier frequency $\omega_c$.


How to recognize hexagonal tiling in boardgame? I would like to recognize the boundaries of a hexagonal tiling in a photograph, like in  image below:

It seems to me, that a standard approach at a square grid is to first detect corners (e.g. canny) and then extract the longest lines via a Hough transform or something similar.

This does not look as optimal solution with hex tiling, because the length of outer lines is shorter and it's hard to separate them from other lines.

Is there any algorithm to adress this problem? It would be particulary nice to have a solution in opencv, but I am also interested in general ideas.

update:

With python and opencv I was able to receive this result:

Here is my code:

import cv2
import numpy as np

imgOrig = "test1";
img = cv2.imread(imgOrig+".jpg");
lap = cv2.Laplacian(img, cv2.IPL_DEPTH_32F, ksize = 3)
imgray = cv2.cvtColor(lap,cv2.COLOR_BGR2GRAY)
ret,thresh = cv2.threshold(imgray,127,255,0)
contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)
size = img.shape
m = np.zeros(size, dtype=np.uint8)
for i, cnt in enumerate(contours):
    if cv2.contourArea(cnt) >= 1:
        color = (255,255,255)
        cv2.drawContours(m, cnt, -1, color, -1)
cv2.imwrite(str(imgOrig)+"contours.jpg", m);

Laplacian of image looks like:

I will try to optimate the parameters of this approach and then try to interpolate the boundaries of the four sections.

1st Approach:

Use the haartraining methods of opencv according to this tutorial http://note.sonots.com/SciSoftware/haartraining.html -- this should give the best results, but I haven't worked with haartraining myself so far...

2nd Approach:

I would suggest to use methods of "markerless tracking" of the individual tiles of the board. You can implement this using OpenCV, too..

Preparation

For this purpose you'll need some photos of each type of tile. Take a picture of all tile types (each one as one picture), with a homogeneous background from top-down-view tile in the middle of the picture.
Then use some feature detector (OpenCV has multiple algorithms for this, but SIFT/SURF are non-free algorithms; I would suggest to use "FAST") to find distinctive points in the images.
Use a feature descriptor to describe the feature found in the image (use e.g. "BRIEF").

Detection

Now you can detect the tiles in an image by applying the same feature detector/descriptor algorithms to this image. When you've acquired the features/descriptors you can apply the FlannBasedMatcher to find the tiles.

Here's a code example / tutorial from OpenCV:
http://docs.opencv.org/doc/tutorials/features2d/feature_homography/feature_homography.html#feature-homography

Notes

The Matcher Method will give you only one match and will possibly have problems if there is more than one tile of that type found on the board. You could work around that problem by masking out only some parts of the input image. I suggest to do this using the pixel coordinates of the detected features. If you - somehow - detect the outline and size of the tiles first, you can roughly estimate the tile positions and size on the picture. Filter your detected feature-list (e.g. only features within x-pixel radius from expected midpoint of the tile) before matching and then use the strongest match. As a result you'll be given the exact position of the tile on the image (including it's orientation). If it's too complicated to detect the map outline, you can let the user "point" at the corner tiles to mark the outline manually...

Alternative Approach

You can also use this method to find just any of the tiles by its outline. Draw a sample "schematic" grayscale picture of a tile (hexagon) without any image on it. Note that the "dark" and "light" regions in this image need to be correct in the schematic, not just some "lines".
You'll probably need to experiment with this. You could try to average multiple photographs of different tiles to generate an "average" image of a tile. Make sure the corners are at the same position (move/scale pictures accordingly) and sharpen the picture when finished (clear corners/edges should be visible) and adjust the contrast a bit, if needed.


Calculate the harmonics in the power supply of $50\textrm{ Hz}$ using FFT I need to calculate the harmonics in the power supply of $50\textrm{ Hz}$.
If suppose I need to calculate up to 20 harmonics.

Then my sampling frequency by Nyquist criteria should be :

$$f_s \ge 2 f_m$$

where $f_m$ will be of the highest harmonic frequency, i.e $f_m = 20 \times 50\textrm{ Hz} = 1000\textrm{ Hz}$.  So $f_s = 2000 \textrm{ samples/sec}$.

If the signal is $50\textrm{ Hz}$ then in this case, $T_{\rm period} = 0.02 \textrm{ sec} = 20\textrm{ msec}$. So in each $20\textrm{ msec}$ period I am taking $0.02 \times 2000 = 40$ samples per period of $50\textrm{ Hz}$.

So in this case do I need to do $40$ point FFT ?
Also if input signal is split into $40$ points then the output of the FFT is also $40$ different values. So this means the output is produced for $40$ harmonics not for $20$ harmonics is it right ?
How this has to be done at software level?

I will take a timer which interrupt at every $1/2000\textrm{ sec} = 0.5\textrm{ msec}$. On every interrupt take the ADC sample. So after $40$ interrupts ($20\textrm{ msec}$ is passed) I will be having buffer of $40$ samples then I will perform the $40$ point FFT algorithm to get the harmonics value, right ?

Then I will clear the buffer & again fill the buffer with new $40$ samples after $40$ interrupts ($20\textrm{ msec}$ is passed) then again I will perform the $40$ point FFT algorithm to get the harmonics value, right ?

As my FFT window is $40$ samples,  do I have to clear the buffer every time after the FFT is performed or I will shift the signal in the buffer to the right by $1$ (i.e. use array of $40$ as ring buffer or FIFO ) ?

Please suggest, is my calculation is right ?

First of all, you should sample at rate greater than 2 times the highest digitized frequency, not equal to.

Second of all, do you have an anti-aliasing filter?


Colorimetry: Gamma/gamut colors Color/color light I've been reading some color basic stuff and I came across with two different terms: gamut and gamma colors. By definition the difference between them may be clear. I found:

Gamma: transfer characteristic of display
Gamut: Range of reproduced color lights

However, these definitions led me to to check these concepts:

Color light: luminance, hue, saturation
Color: Hue and saturation

After checking these last terms everything turned really complex. Could you explain to me what the difference between color and color light is and why the first one does not include the term luminance? How are they connected to gamma and gamut colors?

Thanks in advance :)

As you might know, there are different colour schemes/models available (RGB, HSL, HSV etc) and the important thing to notice is that these schemes are equivalent, which means you can achieve
exactly the same colour manipulations using any of the schemes. Why do we have different schemes then? It's because some things are (much) simpler to do using a scheme designed for a particular purpose.

The hue, saturation and lightness/luminance parameters are a part of the HSL colour scheme. The three parameters together define a colour. It's the luminance parameter that defines the lightness.
This particular scheme  makes it extremely easy to create effects such as colour gradients. The only parameter to change is,of course, the luminance.

As an example, let's say that we want to create a gradient from red to white. In the RGB colour scheme, red is encoded as (255,0,0) and white is (255,255,255). As you can see, both the green and blue components have to be changed
in order to create the gradient.

In contrast, in the HSL scheme red is encoded as (0,255,0) and white is (0,255, 255), which means that we change only the lightness to create our gradient.

To answer your question, all the three parameters constitute a colour whereas only the lighness/luminance defines a colour light.


Fourier Transformed a signal $f(t) = \cos(2\pi*4t)*5\cos(2\pi*0.4t)$ but did not get expected results. What went wrong? I created a signal with two sinusoidal components by specifying

$$
f(t) = \cos(2\pi*4t)*5\cos(2\pi*0.4t)
$$

I expected with I used MatLab's fft function on this that signal that I would see a peak of amplitude $1$ at a position of 4Hz in the spectrum corresponding to the 'fast' signal which has frequency $4$, and a peak of amplitude $5$ at a position of 0.4Hz in the spectrum corresponding to the 'slow' signal which has frequency $0.4$.

In my code below, subplots 3 and 4 are the fft plots, 3 is the unshifted one, and 4 is the shifted one. As you can see from running the code, I do not get the results I expected!

It seems I have peaks at $4+0.4=4.4$ Hz and $4-0.4=3.6$ Hz, and although the scaling is off with the amplitude, both peaks are the same height.

I thought the Fourier Transform would pick up the oscillations at 4Hz and 0.4Hz individually, and display them with their correct (relative) amplitudes? What has gone wrong?

function FFT_4()
clc
clear all
close all

sigFreq = 4;
f = @(t) cos(2*pi*sigFreq*t)*5.*cos(2*pi*0.1*sigFreq*t);

%% Calculation for samFreq > 2*sigFreq: Higher resolution
samFreq = 12;
[t,S,omega,power,fshift,powershift] = CalcFFT(f, samFreq);
plotSignalandFFTSignal(4,t,S,omega,power,fshift,powershift)

end

function [t,S,omega,power,omegashift,powershift] = CalcFFT(f, samFreq)
t = 0:(1/samFreq):(10-1/samFreq); % time vector
S = f(t);
n = length(S);
X = fft(S);
Y = fftshift(X);
omega = (0:n-1)*(samFreq/n);     %frequency range
power = abs(X).^2/n;    %power
omegashift = (-n/2:n/2-1)*(samFreq/n); % zero-centered frequency range
powershift = abs(Y).^2/n;
end

function plotSignalandFFTSignal(figNum,t,S,f,power,fshift,powershift)
figure(figNum)
subplot(2,2,1); plot(t,S); % xlim([0,1])
subplot(2,2,2); plot(power)
subplot(2,2,3); plot(f,power)
subplot(2,2,4); plot(fshift,powershift)
end

Fourier Transform is linear transform. It is mean that Fourier Transforms of sum of functions will be sum of Fourier Transforms of each function. It is not valid for multiply of function.

Your example. There is very well known formula (see 7c in this link):
$$
\cos(\alpha) * cos(\beta) = 1/2 (\cos(\alpha+\beta) + \cos(\alpha-\beta))
$$
use this formula  to transform your function:
$$
f(t) = \cos(2\pi*4t)*5\cos(2\pi*0.4t) = 5/2*(\cos(2\pi*4t - 2\pi*0.4t) + \cos(2\pi*4t - 2\pi*0.4t)) = 5/2*(\cos(2\pi*3.6t) + \cos(2\pi*4.4t))
$$
So Matlab results is correct - you have 2 peaks equal amplitude with 3.6 and 4.4 frequency.


Hf ionospheric channel model for large bandwidth I know some models like the Watterson model for HF ionospheric channels but those are limited to narrow band frequencies.

Are there some models to bandwidth like 12kHz up to 24 kHz?

There is the Watterson channel model, and it can be extended to any bandwidth


Decimating Long Data in Chunks I have a very long data up to 8 millions of samples and want to decimate it to order of 1000s. I'm using following Matlab code to do this:

f = numberOfSamples/N; %N is number of points after decimation
p = factor(f); %
temp = rawData; %with length of numberOfSamples
for i = 1:length(p)
    temp= decimate(temp,p(i),10); %decimate with 10th order iir filter
end

Number of samples before and after decimation is not known, but the ratio, i.e. decimation number, is an integer for sure(that's given.) In other words, those are input to my program.
How can I decimate(including the decimation filter) the data in chunks and have the same result? One chunk computes some part of the decimated signal and the other chunk computes another part, etc.
I hope question is clear, any help is appreciated.

The decimate function in matlab does 2 things. First, the signal is filtered. From your code, and the function doc page, I see that it uses Chebyshev Type I IIR filter of order 10 with normalized cutoff frequency 0.8/p(i) and passband ripple 0.05 dB. Next, the signal is downsampled.
The filter can be created using

[a_lp,b_lp] = cheby1(10,0.05,0.8/p(i));

You should use the initial condition argument in the filter function in order to work in blocks.

Now, the downsampling is a little bit tricky because you have to make sure that the decimation distance is kept between blocks. IMO, the easiest way to do so is to use block size that is multiple of the decimation factor.

Here is something to start with:

N = K*p(i); % K integer
li = 1;
lid = 1;
ri = N;
rid = K;
z = [];
while ri<=numberOfSamples
  [filteredData,z] = filter(a_lp,b_lp,rawData(li:ri,z);
  decimatedData(lid:rid) = filteredData(1:p(i):end);
  li = li + N;
  ri = ri + N;
  lid = lid + K;
  rid = rid + K;
end

(Not sure about the indexes in decimatedData).


How to extract IF (instantaneous frequency) from Discrete Hilbert Transform? The continuous time Hilbert transform is
$$\hat x(t) := x(t) + j\left( p.v. \left\{\frac{1}{t\pi} \ast x \right\} (t)\right)$$

where
$$ \theta(t) = \tan^{-1}\left( \frac{ p.v.  \left\{\frac{1}{t\pi} \ast x \right\} (t)}{x(t)} \right)$$
and
$$
         \omega(t) = \theta'(t)
$$

and since $\theta(t)$ is smooth, you can actually calculate the instantaneous frequency. But in the discrete time domain, the angle function is discrete, so given that in the dt domain, your hilbert transform filter is

$$ h_{HT}[n] = \left\{
      \begin{array}{ll}
       \frac{1}{n\pi}  & \mbox{if}\ \  n \neq 0 \\
        0               & \mbox{otherwise} \\
      \end{array}
     \right.$$

How do you get the approximate frequency? or are you really just getting $\Delta\theta[n]$? Sorry if this is me just getting hung up on something trivial, but most references don't seem to clearly address this.

$$\theta[n] = \tan^{-1}\left(\frac{(h_{HT} \ast x)[n]}{x[n]}\right)$$
$$\omega[n] = \left\{
       \begin{array}{ll}
        \theta[n] \ast \frac{1}{T_s}(\delta[n] - \delta[n-1])  & ? \\
        \theta[n] \ast \frac{1}{2T_s}(\delta[n+1] - \delta[n-1])  & ? \\
        ???
       \end{array}
      \right.
$$

I'm looking for what is typically done. The "Standard" approach.

I have a nice solution for you that I'm always using whenever I need to compute the instantaneous frequency of a discrete analytical signal and which works much better then the common approach using $tan(.)$ where you need to introduce thresholds and fiddle around with adding/subtracting $2\pi$. In addition to that it is even faster.

Given a signal $x(t)$ and the corresponding analytical signal $$\xi(t) = x(t) + i\mathcal{H}x(t)$$ with $\mathcal{H}$ being the hilbert transform operator the instantaneous frequency $\omega(t)$ can be computed as followed. We have $\omega(t)=\dot{\varphi}(t)$ with $\varphi(t)=\Im log\xi(t)$ being the phase of the analytical signal $\xi$ ($\Im$ denotes the imaginary part). Combining these equations we get
$$
\omega(t)=\Im\frac{d}{dt}log\xi(t)=\Im\frac{\dot{\xi}(t)}{\xi(t)}
$$
This can also be seen be as follows:
$$
  \Im\frac{\dot{\xi}}{\xi} = \Im\frac{ \dot{A}e^{i\varphi} + iAe^{i\varphi}\dot{\varphi} }{Ae^{i\varphi}}
                       = \Im\Big[ \frac{\dot{A}}{A} + i\dot{\varphi} \Big]
                       = \dot{\varphi}
$$
Discretizing this equation you can now compute the discrete instantaneous frequency from the discrete analytical signal:
$$ \omega[n]\approx\Im\frac{\xi[n+1]-\xi[n]}{\xi[n]} $$
or second order:
$$ \omega[n]\approx\Im\frac{\xi[n+1]-\xi[n-1]}{2\xi[n]} $$

Of course this is not the exact solution as it is still requires a discrete differentiation of the signal. But this could be done using standard techniques.


How to triangulate 2 sources on a 2d plane formed by 4 sensors? I have 4 sensors capturing signal data. Signal distribution speed is constant and known. I have sensors positions in 3d space.
Also, I have an amplitude/frequency decomposition for data received by each sensor over time (imagine a 3d terrain).
 Sensors capture signal from 2 moving sources.

Here on the image: blue dots are the sensors; red solid Capsules are emmiting objects; transparent red capsules are here to show signal wave propagation pattern; yellow cubes show projections from objects onto a plane formed by sensors; green "points" are desired positions on a plane.

Details:

In real life emitters have pill-shaped form emitting in all directions, yet we can assume that they are points emitting spherical waves.

In my particular case, we discuss audio waves via air to microphones (but I hoped for a general solution).
We do not know the exact signal being transmitted from any of the emitters, yet we do know that there are only two signal sources with not equal positions in space and noise is filtered and extremely low.
We can perform any kinds of transformations on top of received signals and we have exact time-data correlation per each sensor.
Sensors form a square with side of length $l$

So how to triangulate sources positions on a 2d plane formed by 4 sensors - what is the algorithm?

Conceptual partial answer:

You can use 3 mics to find position in 3D space, except you can't tell if it's above the plane or below the plane (unless you already know there's nothing above the plane, of course).

Since your 4 microphones are in the same plane, the 4th microphone is redundant and doesn't help you know whether the source is above or below the plane.  (Though it improves SNR/accuracy of estimate, as would any number of additional mics in the same plane).

There is a delay due to speed of sound from the source ($S$) to each microphone.  The delay difference between each pair of microphones gives you a hyperboloid (green cross-section) that the source must be on, in order for the waves to be delayed by that relative amount (due to extra distance $a$).  Each pair of mics defines a hyperboloid, so by combining 3 pairs of mics, you can find the intersection of the 3 hyperboloids, which will define a unique set of 2 points.  (Symmetrical above and below the plane.)

(4 mics gives you 6 pairs of microphones, so by averaging the intersection point of all 6 hyperboloids you can improve the accuracy of the estimate, but you still don't know if it's above or below the plane.   With a tetrahedron of mics you would only have one intersection point.  You could probably weight the average, too, since some source positions will be more accurately measured by some pairs of mics, etc.)

Also, the sound pressure drops off as 1/r with distance.  The attenuation (amplitude) ratio between each pair of microphones gives you a sphere (yellow cross-section) that the source must be on.  By combining this information from every pair of microphones you can find average intersection of all 6 spheres and improve accuracy from that, too.

The spheres and hyperboloids are usually not tangent at the intersection point, so they each improve accuracy in a particular direction:

To find the relative delay, you can use cross-correlation of the signal from the 2 mics.  There will be a peak at the "lag" at which the two signals line up with each other.

To find the projection of the source location on the plane itself, just throw away the z coordinate.

See https://en.wikipedia.org/wiki/Multilateration#Solution_algorithms for more details on actual implementation.


Determining the type of filter on an image I have the following mask where the origin is marked in bold text

$$ \left[ \frac{-1}{4}~~~{\bf\frac{1}{2}}~~~\frac{-1}{4} \right ] $$

After computing the DFT, the result is:

$$ \frac{1}{2} - \frac{1}{2} \cos(\omega)$$

How can I know what type of filter is it? All-pass, Low-pass, Band-pass, High-pass, Band-cut or All-cut?
To see the type of filter should I just plot it under the main frequency $ [-\pi, \pi] $?

As you said, you can simply plot the frequency response of your filter from $-\pi$ to $+\pi$.

However since this is a filter only has real co-efficients, you can also just plot the frequency response from $0$ to $\pi$, as this will give you all the information you need.

This is because for a real signal, the DFT is conjugate symmetric, meaning that all positive frequencies are mirror images of the negative frequencies. In other words, you do not gain any additional knowledge by looking at both negative and positive frequencies. You just need to look at one of them. (Either negative or positive).

If we do that, we see that your filter looks like so. Here I have shown all positive frequencies, from $0$ to $\pi$. This will give all the information you need.

As you can see, it looks like your filter is a 'graceful high pass filter'. We can clearly see that it passes high frequencies, and clearly see that it nulls low ones, esp DC. In fact, the sum of all your filter co-efficients is 0, meaning we expect it to null out DC completely.

Further judging by context, this appears to be an edge detection filter. Edges are characterized by sharp transitions in space/time, meaning high frequencies in temporal/spatial frequency.

Such a filter would isolate edges.


What is the sparse Fourier transform? MIT has been making a bit of noise lately about a new algorithm that is touted as a faster Fourier transform that works on particular kinds of signals, for instance: "Faster fourier transform named one of world’s most important emerging technologies". The MIT Technology Review magazine says:

  With the new algorithm, called the sparse Fourier transform (SFT), streams of data can be processed 10 to 100 times faster than was possible with the FFT. The speedup can occur because the information we care about most has a great deal of structure: music is not random noise. These meaningful signals typically have only a fraction of the possible values that a signal could take; the technical term for this is that the information is "sparse." Because the SFT algorithm isn't intended to work with all possible streams of data, it can take certain shortcuts not otherwise available. In theory, an algorithm that can handle only sparse signals is much more limited than the FFT. But "sparsity is everywhere," points out coinventor Katabi, a professor of electrical engineering and computer science. "It's in nature; it's in video signals; it's in audio signals."

Could someone here provide a more technical explanation of what the algorithm actually is, and where it might be applicable?

EDIT: Some links:

The paper: "Nearly Optimal Sparse Fourier Transform" (arXiv) by Haitham Hassanieh, Piotr Indyk, Dina Katabi, Eric Price.
Project website - includes sample implementation.

The idea of the algorithm is this: assume you have a length $N$ signal that is sparse in the frequency domain. This means that if you were to calculate its discrete Fourier transform, there would be a small number of outputs $k \ll N$ that are nonzero; the other $N-k$ are negligible. One way of getting at the $k$ outputs that you want is to use the FFT on the entire sequence, then select the $k$ nonzero values.

The sparse Fourier transform algorithm presented here is a technique for calculating those $k$ outputs with lower complexity than the FFT-based method. Essentially, because $N-k$ outputs are zero, you can save some effort by taking shortcuts inside the algorithm to not even generate those result values. While the FFT has a complexity of $O(n \log n)$, the sparse algorithm has a potentially-lower complexity of $O(k \log n)$ for the sparse-spectrum case.

For the more general case, where the spectrum is "kind of sparse" but there are more than $k$ nonzero values (e.g. for a number of tones embedded in noise), they present a variation of the algorithm that estimates the $k$ largest outputs, with a time complexity of $O(k \log n \log \frac{n}{k})$, which could also be less complex than the FFT.

According to one graph of their results (reproduced in the image below), the crossover point for improved performance with respect to FFTW (an optimized FFT library, made by some other guys at MIT) is around the point where only $\frac{1}{2^{11}}$-th to $\frac{1}{2^{10}}$-th of the output transform coefficients are nonzero. Also, in this presentation they indicate that the sparse algorithm provides better performance when $\frac{N}{k} \in [2000, 10^6]$.

These conditions do limit the applicability of the algorithm to cases where you know there are likely to be few significantly-large peaks in a signal's spectrum. One example that they cite on their Web site is that on average, 8-by-8 blocks of pixels often used in image and video compression  are almost 90% sparse in the frequency domain and thus could benefit from an algorithm that exploited that property. That level of sparsity doesn't seem to square with the application space for this particular algorithm, so it may just be an illustrative example.

I need to read through the literature a bit more to get a better feel for how practical such a technique is for use on real-world problems, but for certain classes of applications, it could be a fit.


fft and inverse fft error size and applying impulse response matlab I have a signal called noise. I am using the following command to calculate the mean squared error after fft and inverse fft process:

 sum((noise-real(ifft(fft(noise,2048),num_of_samples))).^2)/num_of_samples

Where num_of_samples=2819519 is the number of samples in the original signal.
I am getting an error of $\sim10^{-8}$. Is that reasonable and due to approximation or does that indicate something is not converted correctly?

My final goal is to multiply the signal in the frequency domain with a $1052\times 4$ given impulse response. This represents a room with 4 microphones.

I don't know your overall goal but if you want to test the MSE associated with FFT / IFFT then you should perform the following

num_of_samples = 2819519 ;
fft_len         = 2^nextpow2(2819519) ;

X = fft(signal, fft_len);
y = real( ifft(X, fft_len) );
mse = sum( (signal - y(1:num_of_samples) ).^2 )/num_of_samples


Power Spectral Density From Variance For an exercise, we made a noise signal by generating a normally-distributed random vector in MATLAB. The mean is defined to be equal to zero and the variance is set to 0.5.

I know that the power spectral density is just the fourier transform of the variance in frequency domain. In this case, does the PSD just equals to 0.5 since the variance is constant with time?

The power spectral density is NOT the Fourier transform of the variance. For wide-sense-stationary (WSS) signals, the power spectral density is the Fourier transform of the Autocorrelation function, and the variance is the autocorrelation at time offset $\tau = 0$ (when the WSS process is zero-mean).
Consider for example a white noise process (which is done if the OP's randomly generated samples were completely uncorrelated from sample to sample), to a tightly filtered process, both can have the same variance (which is the power if WSS and zero-mean), but VERY different spectrums. A white noise spectrum is constant over all frequencies while the filtered process would be limited to a narrow bandwidth. The relationship between the two is given by Parseval's theorem; the sum of the squares of each frequency component will be the same in both cases.
See this answer from Matt L further detailing variance with relation to the signal power: variance in the time domain versus variance in frequency domain
Update:
Dilip stated this all much more elegantly and concisely in the comments, paraphrased here: The area under the power spectral density function of a WSS process equals the value of the autocorrelation function at offset $\tau=0$. For a zero-mean WSS process, this value is the variance. Therefore deducing the power spectral density from knowledge of the variance is equivalent to asking for a curve given that we only know the area, which we know is not possible.


A sequence OFDM modulated with normalization factor The following sequence is used to modulate 52 subcarriers. Twelve out of 52 are non zero. I was wondering

$S_{[–26, 26]} = \sqrt{13/6} \{0, 0, 1+j, 0, 0, 0, –1–j, 0, 0, 0, 1+j, 0, 0, 0, –1–j, 0, 0, 0, –1–j, 0, 0, 0, 1+j, 0, 0, 0, 0,
0, 0, 0, –1–j, 0, 0, 0, –1–j, 0, 0, 0, 1+j, 0, 0, 0, 1+j, 0, 0, 0, 1+j, 0, 0, 0, 1+j, 0,0\}$

The book says that The multiplication by a factor of $\sqrt{13/6}$ is in order to normalize the average power of the resulting OFDM
symbol, which utilizes 12 out of 52 subcarriers.

I tried to compute it..

My guess is that the average energy in such sequence before normalization would be $$12\times \frac{\sqrt{2}}{52}=\frac{3\sqrt{2}}{13}$$

where $\sqrt{2} = \sqrt{1^2+1^2}$ from the complex symbols withing the sequence.

Then why is the normalization factor $\frac{\sqrt{13}}{\sqrt{6}}$?

Thanks

The power of the non-zero elements is $|1+j|^2=2$, so the average power is

$$\overline{P}=\frac{12\cdot 2}{52}=\frac{24}{52}=\frac{6}{13}$$

So in order to get unity average power, you have to normalize by $\sqrt{13/6}$.


Finding the Wiener filter transfer function I am trying to understand the mechanism of finding the transfer function for a Wiener function.

|Rxx(0)   Rxx(-1) ...... Rxx(1-N)| | h(0) |   | Ryx(0) |
|Rxx(1)   Rxx(0)  ...... Rxx(2-N)| | h(1) | = | Ryx(1) |
|......   ......  ......  ...... | | .....|   |....... |
|Rxx(N-1) Rxx(N-2) ...... Rxx(0) | |h(N-1)|   |Ryx(N-1)|

Rxx is the autocorrelation of my input signal h is filter transfer function and Ryx is crosscorrelation of the input and output signal. My question is how to come up with the h. I mean i have two unknowns in this equation cause i dont know h and I dont know the y output signal therefore i cant calculate the crosscorelation Ryx.

Is there something that I need to assume in order to get the h??

Thanks in advance for the help.

What you've listed is a classical set of linear equations known as the Wiener-Hopf Equations. There are many methods to solve this set of equations. These equations generalize the Yule-Walker equations for Autoregressive (AR) modeling. In fact, if the Wiener-Hopf equations are used to solve for a linear predictor of your signal, the equations are identical to the Yule-Walker equations.

So in reality, there are many ways to solve these equations depending on how you're modeling your signal. This is where the concepts of Wiener filtering, Wiener smoothing, Wiener deconvolution arise. In each of these, they can be segmented further into causal/anti-causal realizations, FIR/IIR, etc. I recommend searching for the Wiener-Hopf equations and understand what they mean, then solving them comes secondary.

As an example, consider if you assume that you have a signal with additive noise that is uncorrelated with the signal. Then the cross-correlation terms have no dependence on the signal. Then you need to model your noise and plug in the resulting autocorrelation sequence on the right hand side of the equation.

Cheers


Is there a way to obtain the impulse response of a discrete system by just knowing it's response to the discrete unit step function? In continuous time it was possible;

$$
u(t){\longrightarrow} \boxed{\quad\textrm{system}\quad} {\longrightarrow} y(t)\implies
\delta(t)=\frac{du(t)}{dt}{\longrightarrow}\boxed{\quad\textrm{system}\quad}{\longrightarrow} \frac{dy(t)}{dt}=h(t)
$$

Does the same apply for discrete time system i.e.
$$
\delta[t]=\frac{du[t]}{dt} \quad\textrm{where:}\begin{cases} \delta[t] &\textrm{is the discrete time delta}\\
u[t] & \textrm{is the discrete time unit step function}\end{cases}
$$

Is there a way to obtain the impulse response of a discrete system by just knowing the response of the discrete unit step?

A simpler version of Phonon's answer is as follows.

Suppose that $y$ denotes the response of the system to the unit step function.
Then, as discussed in this answer,
in general, $y$ is the sum of scaled and time-delayed copies of the impulse
response, and in this particular case, no scaling is required; only time
delays.  Thus,
$$\begin{align}
y[0] &= h[0]\\
y[1] &= h[1] + h[0]\\
y[2] &= h[2] + h[1] + h[0]\\
y[3] &= h[3] + h[2] + h[1] + h[0]\\
\vdots~ &= ~\vdots
\end{align}$$
where each column on the right is a (unscaled and) time-delayed
impulse response.
Thus, we easily get that
$$\begin{align}
h[0] &= y[0]\\
h[1] &= y[1]-y[0]\\
h[2] &= y[2]-y[1]\\
\vdots~ &= ~\vdots\\
h[n] &~= y[n] - y[n-1]\\
\vdots~ &= ~\vdots
\end{align}$$
with nary a mention of filters, inverses, convolutions, integration, operators
and the like, just simple consequences of
the definition of linear time-invariant system.


Calculating distance between two blobs in a binary image in microns I have a monochrome CCD camera with resolution 960 X 1280 and can give up-to 15 fps
I am working on python opencv.

I want to track the distance between two rectangular blobs in binary images in real time (1 value per 2 or more secs)

My field of view is 40 mm and i want to obtain the resolution of upto 1 micron with +-1% accuracy.

Can this be achieved using binary image or is there any other way to achieve this ?

all suggestions are welcome

This is relatively straight-forward with OpenCV and numpy.

Step One:
Using OpenCV's Blob Detector you can do something like this:

import cv2

im = cv2.imread("blobs.png", cv2.IMREAD_GRAYSCALE)
params = cv2.SimpleBlobDetector_Params()
params.blobColor = 255
detector = cv2.SimpleBlobDetector_create(params)

keypoints = detector.detect(im)

kpt_image = cv2.drawKeypoints(im, keypoints, numpy.array([]),(0,0,255),cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

cv2.imshow("Keypoints", kpt_image)
cv2.waitKey(0)
cv2.destroyAllWindows()

Step Two:
Write a function to calculate the distance between two keypoints:

import numpy

def distance(kpt1, kpt2):
    #create numpy array with keypoint positions
    arr = numpy.array([kpt1.pt, kpt2.pt])
    #scale array to mm
    arr = arr*40/1280
    #return distance, calculted by pythagoras
    return numpy.sqrt(numpy.sum((arr[0]-arr[1])**2))

Putting it all together:

import cv2
import numpy

def distance(kpt1, kpt2):
    #create numpy array with keypoint positions
    arr = numpy.array([kpt1.pt, kpt2.pt])
    #scale array to mm
    arr = arr*40/1280
    #return distance, calculted by pythagoras
    return numpy.sqrt(numpy.sum((arr[0]-arr[1])**2))

im = cv2.imread("blobs.png", cv2.IMREAD_GRAYSCALE)
params = cv2.SimpleBlobDetector_Params()
params.blobColor = 255
detector = cv2.SimpleBlobDetector_create(params)

keypoints = detector.detect(im)

kpt_image = cv2.drawKeypoints(im, keypoints, numpy.array([]),(0,0,255),cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

for i,keypoint in enumerate(keypoints[1:]):
    print("Distance: {0:6.3f}mm".format(distance(keypoints[0], keypoint)))

cv2.imshow("Keypoints", kpt_image)
cv2.waitKey(0)
cv2.destroyAllWindows()

I assume that when you say 'real time' you may want to look at annotating the image with the distance.  But I'll leave you to work with that yourself.  I've not actually done an error analysis, but you may need to tweak the Params above to increase accuracy.

The output image:


Minimum Number of Errors required for BER Curves in Matlab Simulation what should be minimum Number of Errors required for getting BER Curves in Matlab Simulation for different number of transmitted bits ?

I have done a fair number of BER simulations and have not found a definitive answer on this. It depends on the error behavior of your channel, Gaussian or some other pdf. My process for AWGN was as follows. First I determine the approximate BER I am expecting. Say that it is 10-3. The inverse of that is 1000 bits. So I would run at least 10x that to start, or a length of 10,000 bits as a minimum. For coded signals, or those with fading, I might use a larger multiplier, such as 50 to 100 times. By watching the error pattern we can often tell if steady state has been achieved. For coded channels with BERs of 10-8 etc., the run lengths quickly become an issue.

This applies to Monte Carlo simulations. You can do semi-analytic methods that do not require such large runs. However semi-analytic methods are generally not that good for non-AWGN channels. Importance sampling is one such semi-analytic method which is used with coded signals. All of these are estimates.


How to calculate the Gain (in dB) given ω₁ and sampling time, T, of a discrete transfer function? Having the Transfer Function of a discrete system as such:

$$H(z) = (1.011)\frac{z-0.9756}{z-0.9975}$$

How would you find the gain at ω₁=3 rad/s at a sampling time of T=0.25 seconds.

The answer to the solution is: 0.1013 dB

I'm having problems with this. I am aware that $z=e^{sT}=e^{jωT}$. I substitute the values of ω₁ and T, and it always results in a complex error, despite attempting it multiple times.

Could someone help how you would get the final solution, with steps provided?

Kind regards,
Alan.

First of all, convert the Z-transform $H(z)$ into the Frequency response $H(e^{j\omega})$ by replacing $z$ with $e^{j\omega}$.

$$ \begin{align}
H(z) &= (1.011) ~\frac{z-0.9756}{z-0.9975} \\ \\
&= (1.011)~\frac{1-0.9756 ~z^{-1}}{1-0.9975~z^{-1}} \\ \\H(e^{j\omega}) &= (1.011)~\frac{1-0.9756 ~e^{-j\omega}}{1-0.9975 ~e^{-j\omega}}\\
\end{align}
$$

Then find the discrete-time radian frequency $\omega$ from the given continous-time frequency $\Omega = 3 ~\text{rad\s}$ by the  sampling relation (assuming no aliasing):

$$ \Omega = 3 ~ \text{rad/s} \implies \omega = \Omega \times T_s = 3 \times 0.25 = \frac{3}{4} ~\text{rad/sample}$$

Then find the gain of the LTI system at the given frequency $\omega = 3/4$ :

$$
H(e^{j3/4}) = (1.011)\frac{1-0.9756 ~e^{-j3/4}}{1-0.9975 ~e^{-j3/4}}
$$

Gain is the magnitude: $|H(e^{j3/4})| = 1.0004$. And the gain in dB is given by:

$$ G_{db} ~= 20 \log_{10}( |H(e^{j3/4})| ) ~= 0.0035 ~~\text{dB} $$


What happens when we Fast Forward the Video I've been involved in Image Processing kind of project for the last couple of month, before which i had no experience in Image processing.
While filtering the image today i had a thought what exactly happens when we fast forward any movie or a video clip ?

For example,
If I have a video whose Frame rate is 30 , this means that 30 frames will be process in a second.
and each frame will take 1/30 = 0.03sec to process , If we increase the frame rate or FPS of the video we will get more smooth view of it.

So what exactly happens when we fast forward the video, what i think is (not sure) the audio and the video of any clip must be synchronized so that we can control both accordingly.

So when we forward the video does it skip certain Frames ? if yes, then what about the Audio ?
What exactly happens ?

What happens to the video:

In the most simple implementation, which is suitable for large integer fast forward rates, only 1 frame out of N, where N is the fast forward rate is played. For example, with a normal play rate at 30 FPS, frames 0 to 29 are played in sequence over a duration of one second. In fast forward with a 3x rate, frames 0, 3, 6 ... 87 are played in one second; and frames 1, 2, 4, 5... 88, 89 are discarded. This can be quite computationally expensive as the video decoder has to provide 3 times as many frames per second as in normal play (though the display continues receiving 30 frames per second).
When the hardware permits it and/or when the fast forward rate is small, the frame rate is simply increased. For example, in fast forward with a 1.1x rate, frames 0 .. 32 are played in one second.
There are applications when a video stream has to be sped up/down by a small factor (say 5%) while maintaining a target frame rate (for example when a broadcaster wants to speed up a show to make room for twenty extra seconds of ads). In this case, interpolation/resampling methods have to be used - this is the same kind of techniques as those used for telecine conversion.

What happens to the audio:

In some implementations, the audio is read at a higher sample rate, and decimated to the target sample rate. For example let us consider the example where the original audio is at 44.1kHz, and the fast forward rate is 3x. In 1s, you need to play 3s of the original audio. So you read 44100 x 3 samples from the original audio stream, low-pass filter it with a filter which has a 1/3 normalized cutoff frequency; and keep only every third sample. This gives you 44100 samples of sped-up audio, which are played in a time span of 1s at 44.1kHz. This has the unwanted effect of shifting up the pitch of all sounds, giving voices a "chipmunk effect". For small shift up/rates the effect is acceptable to the common public (though musically trained people might notice a well-known song has been shifted in pitch and is not in its original tonality).
In some better implementations (eg Apple Quicktime player), a time-stretching algorithm is used to multiply the speed of the input audio without altering the pitch. For playback / scrubbing, the sound quality rarely matters and the time-stretching algorithms are rather crude (you can have a look at the source code of soundtouch).


Phase shifting a noisy signal I have a signal of the form $s(t)=A(t) \sum \cos(\omega_i(t)t +\phi_0) + n(t)$, where $n$ is gaussian noise.

I can only read the signal+noise and thus can not separate them.

I want to phase shift the signal to $A(t) \sum\cos(\omega_i(t)t)+n'(t)$ and I am at a loss on how to do this. During lectures / courses I've always done phase shifts simply by multiplying with $e^{i\phi_0}$.

My signal is of the form $s(t)=A(t) \sum \cos(\omega_i(t)t +\phi_0) + n(t)$, as opposed to $s(t)=A(t) \sum e^{\omega_i(t)t +\phi_0} + n(t)$, so I can not simply multiply with $e^{i\phi_0}$

Is there any way to do this? I'm asking because I am interested in cross-correlating two signals.

With a real cosine a phase shift is equivalent to a time delay.  If you delay the signal by $\phi_0$ then you will get the result you want.

The number of samples that you will need to delay by will likely not be an integer number.  You can delay by a non-integer number of samples via a fractional delay filter.


Notch filters: why is $H(e^{j0}) = H(e^{j\pi})$ required? I'm studying for the final exams and many example exercises require to design fir notch filters of length 5, given some passband and some stopband frequencies. You've got to find 3 coefficients α0, α1, α2. You need 3 equations.

R(passband ω) = 1
R(stopband ω) = 0

I completely understand why these are here. The final degree of freedom is determined by the equation $H(e^{j0}) = H(e^{j\pi})$, because such thing is true for notch filters. I struggle to understand why this is true.

In general, there is no such requirement for notch filters that $H(e^{j0})=H(e^{j\pi})$ must be satisfied. You could definitely have a notch filter with $H(e^{j0})\neq H(e^{j\pi})$. Having the same gain at DC and at Nyquist is just a practical definition, and if you have a sufficient number of degrees of freedom (i.e., filter coefficients) you might as well make it a requirement. Ideally you want $|H(e^{j\omega})|=1$ for all frequencies except for the notch frequency $\omega_0$, for which you require $H(e^{j\omega_0})=0$.

As a simple example, take the most primitive notch filter, namely a FIR filter of length $3$. Its frequency response must be

$$H(e^{j\omega})=c\cdot \left(1-2z^{-1}\cos(\omega_0)+z^{-2}\right)\tag{1}$$

where $\omega_0$ is the notch frequency, and $c$ is a scaling constant. All degrees of freedom are taken by the requirement of the notch at $\omega_0$, so you don't get equal gains at DC and at Nyquist. Yet, it is a notch filter.


Is there any alternative characterization of sparsity of a signal in compressed sensing The starting assumption for compressed sensing (CS) is that the underlying signal is sparse in some basis, e.g., there are a maximum of   non-zero Fourier-coefficients for an  $s$-sparse signal. And real life experiences do show that the signals under consideration are often sparse.

The question is - given a signal, before sending out the compressively-sampled bits to the receiver and let her recover to the best of her abilities, is there a way to tell what its sparsity is, and if it is a suitable candidate for compressed sensing in the first place?

Alternatively, is there any additional/alternative characterization of sparsity that can tell us quickly whether CS will be useful or not. One can trivially see that the sender could do exactly what the receiver will do with some randomly chosen set of measurements, and then try to figure out the answer. But is there any alternate way to resolve this question ?

My suspicion is that something like this must have been studied, but I couldn't find a good pointer.

Note : I had posted this question in Mathoverflow, a few weeks back, but didn't get any answer. Hence the cross-post.

Indeed, there are ways in which sparsity, or information content, may be estimated
at the acquisition device. The details, practicality, and actual usefulness of doing so
is debatable and heavily dependent upon the context in which it is applied. In the case
of imaging, one could determine areas of an image which are more or less compressibile
in a predetermined basis. For example, see "Saliency-Based Compressive Sampling for Image Signals" by Yu et al. In this case, the additional complexity requirements placed on the acquisition device provide marginal gains.

With regards to your questions about making determinations as to the usefulness of Compressed Sensing on a given signal at the time of acquisition: If the signal in question adheres to
any kind of model known a priori, Compressed Sensing is possible. Accurate recovery is
simply dependent upon ratio between the number of measurements taken and the degree to which
the sampled signal adheres to your model. If it is a bad model, you won't get past the phase
transition. If it is a good model, then you will be able to calculate an accurate
reconstruction of the original signal. Additionally, Compressed Sensing measurements are,
in general, future proofed. If you have a given number of measurements for a signal which
are insufficient in number of accurately recover the original signal using the model you
have today, then it is still possible to devise a better model tomorrow for which these
measurements are sufficient for accurate recovery.

Additional Note (edit):
The acquisition approach mentioned in your question sounded quite near to adaptive Compressed Sensing, so I thought the following might be of interest for the readers of this question. Recent results by Arias-Castro, Candes, and Davenport have shown that adaptive measurement strategies cannot, in theory, offer any significant gains over non-adaptive (i.e. blind) Compressed Sensing. I refer readers to their work, "On the Fundamental Limits of Adaptive Sensing" which should be appearing in the ITIT soon.


Profile Matching in a Point Cloud A point cloud is generated using uniform random function for (x,y,z). As shown on the following figure, a flat intersecting plane (profile) is being investigated that matches as the best (even if not the exact one) a target profile i.e., given at the left-bottom corner. So the question is:

  1- How to find such a match of given target 2D point map through point cloud considering the following notes/conditions?
  2- What are then the coordinates/orientations/degree of similarity etc?

Note 1: The profile of interest could be anywhere with any rotation along axes and also could be in different shape e.g., a triangle, rectangle, quadrangle etc depending on its location and orientation. At the following demonstration only a simple rectangle one is shown.

Note 2: A tolerance value could be considered as the distance of the points from the profile. To demonstrate this for the following figure suppose a tolerance of 0.01 times the smallest dimension (~1) so tol=0.01. So if we remove the rest and project all remaining points on the plane of the profile being investigated then we will be able to check its similarity with the target profile.

Note 3: A related topic could be found at Point Pattern Recognition.

This is always going to require a lot of computation, especially if you want to process as many as 2000 points.  I'm sure there are highly-optimized solutions for this kind of pattern-matching already, but you have to figure out what it's called in order to find them.

Since you're talking about a point cloud (sparse data) instead of an image, my cross-correlation method doesn't really apply (and would be even worse computationally).  Something like RANSAC probably finds a match quickly, but I don't know much about it.

My attempt at a solution:

Assumptions:

You want to find the best match, not just a loose or "probably correct" match
Match will have some small amount of error due to noise in measurement or computation
Source points are coplanar
All source points must exist in target (= any unmatched point is a mismatch for the entire profile)

So you should be able to take a lot of shortcuts by disqualifying things and decrease computation time.  In short:

pick three points from source
search through target points, finding sets of 3 points with the same shape
when a match of 3 points is found, check all the other points in the plane they define to see if they're a close match
if more than one match of all points is found, choose the one with the smallest sum of 3D distances error

More detailed:

pick a point from the source for testing s1 = (x1, y1)
Find nearest point in source s2 = (x2, y2)
d12 = (x1-x2)^2 + (y1-y2)^2
Find second nearest point in source s3 = (x3, y3)
d13 = (x1-x3)^2 + (y1-y3)^2
d23 = (x2-x3)^2 + (y2-y3)^2

for all (x,y,z) test points t1 in target:
    # imagine s1 and t1 are coincident
    for all other points t2 in target:
        if distance from test point > d12:
            break out of loop and try another t2 point
        if distance ≈ d12:
            # imagine source is now rotated so that s1 and s2 are collinear with t1 and t2
            for all other points t3 in target:
                if distance from t1 > d13 or from t2 > d23:
                    break and try another t3
                if distance from t1 ≈ d13 and from t2 ≈ d23:
                    # Now you've found matching triangles in source and target
                    # align source so that s1, s2, s3 are coplanar with t1, t2, t3
                    project all source points onto this target plane
                    for all other points in source:
                        find nearest point in target
                        measure distance from source point to target point
                        if it's not within a threshold:
                            break and try a new t3
                        else:
                            sum errors of all matched points for this configuration (defined by t1, t2, t3)

Whichever configuration has the least squared error for all other points is the best match

Since we're working with 3 nearest neighbor test points, matching target points can be simplified by checking if they are within some radius.  If searching for a radius of 1 from (0, 0), for instance, we can disqualify (2, 0) based on x1 - x2, without calculating the actual Euclidean distance, to speed it up a bit.  This assumes that subtraction is faster than multiplication. There are optimized searches based on more arbitrary fixed radius, too.

function is_closer_than(x1, y1, z1, x2, y2, z2, distance):
    if abs(x1 - x2) or abs(y1 - y2) or abs(z1 - z2) > distance:
        return False
    return (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 > distance^2 # sqrt is slow

3D Euclidean distance is $d = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2+(z_1 - z_2)^2}$
but square roots are slow. Since we're just comparing relative distance, I think we could drop the sqrt and consider "distance" to be sums of squares.

The minimum calculation time would be if no 2-point matches are found.  If there are 2000 points in the target, this would be 2000 * 2000 distance calculations, though many would be disqualified by a subtraction, and the results of previous calculations could be stored so you only have to do $2000 \choose 2$ = 1,999,000.

Actually, since you will need to calculate all of these anyway, whether you find matches or not, and since you only care about nearest neighbors for this step, if you have the memory it's probably better to pre-compute these values using an optimized algorithm.  Something like a Delaunay or Pitteway triangulation, where every point in the target is connected to its nearest neighbors.  Store those in a table, then look them up for each point when trying to fit the source triangle to one of the target triangles.

There's a lot of calculations involved, but it should be relatively fast since it's only operating on the data, which is sparse, instead of multiplying lots of meaningless zeros together like cross-correlation of volumetric data would involve.  This same idea would work for the 2D case if you found the centers of the dots first and stored them as a set of coordinates.