Is there a way to reduce the covariance matrix of several source signals to the dominant source signal? The problem I have can be seen in the context of DoA estimation or blind source signal separation and similar fields, where several source signals are observed by several antennas (or by an antenna array, respectively).

The basic model is something like
$$\mathbf{y}(t) = \sum_{i=1}^{n_s} s_i \mathbf{a}\left(\theta _i\right) + \mathbf{n}\left(t\right),$$
where $s_i$ describes the source signals, $\mathbf{a}\left(\theta _i\right)$ the mixing vector and $\mathbf{n}\left(t\right)$ the additive noise term. Note that neither of these terms should be considered known.

From this, I have available the sample covariance matrix
$$\mathbf{\hat{C}} = E\{\mathbf{y}(t)\mathbf{y}^H(t)\} = \frac{1}{T} \sum_{t=1}^T\mathbf{y}(t)\mathbf{y}^H(t)$$
with dimensions $N \times N $, if $N$ is the number of antennas, $n_s$ the number of source signals and $T$ the number of samples.

Now, what I want to do is the following:
I don't want to estimate all DoAs nor all source signals. I am looking for a method that is able to reduce the covariance matrix containing information about $n_s$ source signals to a covariance matrix containing just information about the dominant source. Actually, it is not even so important that it is the dominant source, the main goal is the reduction to just one signal.

Is there any method to achieve this kind of cleaning, of signal elimination, of source reduction, of clutter elimination or whatever you want to call it?

EDIT: I think I have not made clear a very important side fact: I am not particularly interested in the amplitudes or powers $\mathbf{s}$ themselves, but more in their geometry/directions $\mathbf{\theta}$.

What you want to do is dimension reduction. The most basic, yet very powerful and commonly used, technique to do it is principal component analysis (PCA). PCA operates on the covariance matrix. You can look it up on Wikipedia for a throughout tutorial.

PCA decomposes your data $y$, having dimension $d$, into $d$ principal components. The first principal component explains the most variance in your data, the next component the second most variance, and so on. Throw off the last $k$ principal components $-$ and ta-daa $-$ you reduce dimension. To answer your question, take $k = d - 1$ and then you get the most dominant component (in terms of explained variance).

(Tip: MATLAB and similar software usually come with an easy-to-use function to perform PCA).


Implementation of CIC filter using floating point I am trying to implement a CIC filter and was reading through the literature and found that CIC filter cannot be implemented using floating point. I am unable to understand why that would happen. The reasoning that is given everywhere seems to be that the integrators are not stable when using floating point but I am not able to understand why that would be. Can someone explain how that would happen?
I also found that using polyphase decomposition would allow for floating point implementation but I am not sure how that would happen since the transfer function is still the same. I am not able to find a satisfactory explanation for this. Is this the only way to do floating point implementation?
Also I read that CIC filter is good if the decimation factor is high, but I am looking for a decimation factor of about 5 to 20. Why wont using a CIC filter be a good approach in this case.

Heh.  I've just worked this out, and I think you can implement a CIC using floating point -- but it wouldn't be perfect, and in most processing environments it would be a waste of resources.
The basic CIC stage is an integrator implemented modulo a maximum number -- in an integer implementation, this would be $x_{max} = 2^n$ where $n$ is the word length, in a non-integer fixed-point implementation it would be the maximum span of values that the data type could hold.
$$x_k = \left(x_{k-1} + u_k \right) \mod x_{max} \\
y_k = \left(x_k - x_{k - N} \right) \mod x_{max}$$
where "$y = x \mod x_{max}$" maps $x$ onto the interval $y \in \left[ -\frac{x_{max}}{2}, \frac{x_{max}}{2} \right )$.
You choose $x_{max}$ so that over $N$ samples, there's no chance that the value of $x_k$ can roll over more than once.
The thing that makes the CIC less than practical in a floating point environment is that you have to do the modulo operation explicitly, either by actually doing a floating point modulo, or by detecting that $x_k$ or $y_k$ is outside of $\left[ -\frac{x_{max}}{2}, \frac{x_{max}}{2} \right )$ and adjusting them by adding or subtracting $x_{max}$.
The thing that makes the CIC less than exact in a floating point operation is that floating point arithmetic truncates more or less depending on the exponent.
The reason that the CIC is so attractive for fast, computationally efficient DSP is because with fixed point (integer or not) math, is because with typical adder hardware and 2's part arithmetic, the modulo operation comes for free -- it's just the fixed-word-width overflow that's usually the bane of fixed-point arithmetic.  So you pay nothing -- in a sense, less than nothing* -- for the modulo.
The reason that CIC is exact is because if all the arguments are fixed point with the same radix (i.e., integers), then all of the operations are bit-exact.  This can't be said for floating point.
So -- in a floating point environment you can implement a CIC filter.  But, you'll be larding on so much extra computational resources to do the two modulo operations, that you may as well do something else, like half-band filters.  The only place where it would make sense to use a CIC in floating point, in my opinion, would be if you're working on a test bench and you want to simulate a system that you're going to implement later.
* because you don't have to worry about detecting it, or preventing it.


Detecting a fixed template image out of a semi-constant frame of video There are a number of videos that I'm looking to process of different video games to detect various "states" in them.

The first game that I am tackling is any edition of Super Street Fighter 4.

In it, I'd like to detect when the "vs" character screen comes up.  Here's an example of one frame of the video:

(taken from ~10s mark of this video)

If I could detect the "vs", then I'd be able to detect that frame of video is indeed the "vs" screen, which would allow me to look for other information (for now, let's say I'll use it to detect the timestamp in the video where the match is about to start).

That said, here is what can be assumed about the frames from the videos that I will be processing (this is not the only video, there are thousands, if not tens or hundreds of thousands of videos, but the issue of scale in processing that many videos is a completely different domain):

I'd prefer (but it's not necessary) to process the lowest resolution image possible with reliable results (lower resolutions = faster processing time). The image above is 480 x 270 pixels (taken from a YouTube video with a fmt 18) but they may come in different sizes (I've gotten YouTube videos with fmt 18 but with dimensions 640 x 360 pixels).
Most videos will be direct-feed
Most videos will be 16:9 aspect ratio
The reddish background will be animated, but generally be within that orange-red color (it's flames)
Sometimes there will be a badge fading in and out over the lower part of the "vs" to indicate a version (that will be important, but not right now), which might obfuscate the "vs", like so:

(taken from ~3s mark from this video; also note that the above is a resolution of 640 x 360 pixels)

The size and position of the "vs" is going to be roughly the same (I haven't verified this yet but I know it doesn't move) in proportion to other direct-feed videos
The characters will be chosen from a pool of more than 30 on each side (in other words, those areas of the frame will vary)
The videos will generally be anywhere from two to four minutes long, with somewhere between 4,000 and 6,00 frames.  However, there might be longer videos (maybe a two hours) which have various other games and live action cut in.  These videos are not as important, but if a solution tells me where a certain game pops up in the larger overall video, great
The native resolution of the captures is 720p, so a baseline image of the "vs" can be taken at what would be considered a "native" size.

Ultimately, I'm looking to code this pipeline in .NET, but that's not super important, the proof-of-concept is more important here as well as an understanding of the techniques involved so that I can translate and optimize it for .NET as well as for other videos of other games in the same genre (if I can pick out the significant discriminators, and videos of say, Ultimate Marvel vs. Capcom 3, Street Fighter x Tekken, BlazBlue: Continuum Shift, etc.).

I'm also dipping my toes in Mathematica and have home version 8.0, so a proof-of-concepts in that environment is more than welcome as well.

If the "VS" is pretty much the same (save for some badge overlays as in the second example), you can use straightforward cross-correlation to detect the presence of the template in your video frame. I answered a similar question on doing this in MATLAB on Stack Overflow. You can use something like the "magic wand" tool in Photoshop to select the "VS" from the frame to create a template. I've done so and binarized the image to obtain this template.

Looking at the different color channels (RGB) in your two images, the red channel appears to be the best for detecting your template.

You can now cross-correlate the red channel with your binarized template and you should get a peak at the location of the template. I choose to threshold and binarize the red template too, although you can detect it without doing so. I prefer to use a distance function rather than raw cross-correlation values, as it tends to be a bit more robust against false positives. I don't know C#/.NET, but here's an outline of the approach in Mathematica:

image = Import["http://i.stack.imgur.com/7RwAh.png"];
ImageCorrelate[ Binarize[ColorSeparate[image][[1]], 0.1], vsTemplate,
   NormalizedSquaredEuclideanDistance] // Binarize[#, 0.2] & // ColorNegate

which gives you the following. The white dot marks the region with the minimum distance in each of the two images

You can then use the above in your next step as desired. Note that typically, cross-correlation will result in an overhang. In other words (using a 1D example) if you're cross-correlating an $N$ point signal with an $M$ point one, you'll get a result that's $N+M-1$ point long. Mathematica's implementation takes care of the overhang for you. However, I don't know what C# does, and you might want to keep this in mind (MATLAB doesn't do it, and I had to account for this in my linked answer above).

You can also build upon this and implement a more robust thresholding criterion of your own. For now, I shall just highlight the detection for the benefit of others:

You can generate the above with a combined function:

detectVS[i_Image] :=
 Module[{mask =
    ImageCorrelate[ Binarize[ColorSeparate[i][[1]], 0.1], vsTemplate,
       NormalizedSquaredEuclideanDistance] ~Binarize~ 0.2 //
     ColorNegate},

  ColorConvert[i, "Grayscale"]~ImageAdd~
   ImageMultiply[i, Image[mask]~Dilation~ DiskMatrix@100]
  ]

There is a lot of potential for improvement here. I'm an armchair hobbyist at image processing, so I don't know what the fastest algorithms are. However, there are a few things you could look into:

If the VS is is roughly the same location in every video, you needn't cross-correlate using the entire image – you can just select a box in the middle and work with that.
This might be an expensive operation to do for each and every frame. However, looking at your video, you have about a little over 4 seconds worth of frames where you have the VS displayed and the character names. So I would suggest that you analyze a frame every second or at most every 2 seconds, thereby guaranteeing that you will land on one with a VS on it. Once you detect VS, you can then start processing every successive frame to do the next part of your processing.
This process should be, to a reasonable extent, robust towards change in size i.e., you could do cross-correlations on small images, but you'll need a suitable template to match. If you know that your images are going to be in certain set/standard sizes, then you could create templates for each of them and select the appropriate template depending on the image size.
The thresholds I chose were by trial and error, but they do seem to work for the two images above and from the other related youtube videos, they probably will work for most of them. A more specialized approach would involve dividing it up into blocks and looking at the histogram to infer if it belongs to VS or not – perhaps a Bayesian classifier. However be absolutely sure that you need to do this before embarking on it. It seems to me that it's simple enough that you won't need it.


Recovering data from mystery signal I was scanning the airwaves using my RTL-SDR and found a signal being broadcast locally that I was able to receive. When I configure the SDR Console V3 software to NFM, I get the following audio. (Wav file)

So I started investigating and here are my findings:

What I know:

According to my local FCC equivalent the license documents list it as coming from an electrical switch yard.

There is are two different frequencies for Tx and Rx.

The signal runs continuously.

What I assume:

That this is a 2-fsk/BFSK signal, due to the two tones it uses.

That it is broadcast in NFM. That mode does sound the best however, and lines up with licensing documents.

That it is actually data, not just some tone control system or random noise.

What I have done

In SDR Console V3 I created an IQ recording. (Large File)

I then created a flow graph in GNURadio to record a file using a sink file. (Large File)

I created a GNURadio flow graph using that file as an input to low pass filters and a FM Demod.

The Demodulated file was then opened in Inspectrum which produced the following image.

What I think I need to do next

I think I need to recover the clock of the signal somehow. I have had a look for a preamble but can't see anything that might match.

I thought it might be Manchester encoded but I can see long runs of 0's or 1's so I don't think it is.

I thought it might be clocked to the frequency of the local power supply given that it is in a switch yard, but it doesn't line up completely with rise or falls of the signal.

Am I on the right path with this?

Is there a way to determine the start of the message? The receiving radio has to be able to lock on to the start of the message somehow, right?

Does this get you anywhere? Apologies for my lacking coding vocabulary, feel free to edit.

Figure 1. Adobe Audition spectrogram of the lowest band of 06-Sep-2018 18h48m47s - 460.519.5 MHz, NFM.wav using Kaiser (180dB) window, 2048 bands, 100 % window width, 80 dB range. A manually adjusted grid is superimposed over the spectrogram.

It's repeating a 16-bit pattern ???? ???? ?111 1110 (spaces added for readability). That would be 9 data bits, and a synchronization sequence of 6 up bits and 1 down bit. The data bits should not contain the string 1111110 or else synchronization may fail. You can synchronize the bit clock to the first up→down transition after a long enough run of the up frequency.


Autocorrelation of LTI system output I'm having trouble showing the following relation:

The autocorrelation of a LTI system with impulse response $h[n]$, input $x[n]$ and output $y[n]=h[n]*x[n]$ is given as:

$r_{yy} = r_{hh}*r_{xx}$

My approach was the following:

$r_{yy} = E(y[n]y[n+m]) = E((h[n]*x[n])\cdot (h[n+m]*x[n+m])) = E((h[n]h[n+m])*(x[n]x[n+m])*(h[n]x[n+m])*(x[n]*h[n+m])) $

but I don't know how to simplify this further...

Actually, we have two kinds of autocorrelation functions. One is defined for stochastic signals and the other for deterministic ones.

If $x(t)$ is a stochastic signal, then its autocorrelation function will be $$R_x(t+\tau,t)=\Bbb E\{x(t+\tau)x^*(t)\}$$where $\Bbb E(\cdot)$ denotes the mathematical expectation.

If x(t) is a deterministic signal, then its autocorrelation function will be $$R_x(t)=x(t)*x^*(-t)=\int_{-\infty}^\infty x(\tau)x^*(t+\tau)d\tau$$
where $\cdot *\cdot$ denotes the convolution operation.

Based on the previous category, we have four possible cases:

  (1) Both $x(t)$ and $h(t)$ are deterministic.

In this case, $y(t)$ is deterministic and we obtain $$R_y(t){=\int_{-\infty}^\infty y(\tau)y^*(t+\tau)d\tau\\=\int_{-\infty}^\infty \int_{-\infty}^\infty\int_{-\infty}^\infty x(\tau_1)h(\tau-\tau_1)x^*(\tau_2)h^*(t+\tau-\tau_2)d\tau_1d\tau_2d\tau\\=\int_{-\infty}^\infty \int_{-\infty}^\infty x(\tau_1)x^*(\tau_2)\int_{-\infty}^\infty h(\tau-\tau_1)h^*(t+\tau-\tau_2)d\tau d\tau_1d\tau_2\\=\int_{-\infty}^\infty \int_{-\infty}^\infty x(\tau_1)x^*(\tau_2)\int_{-\infty}^\infty h(\tau)h^*(t+\tau+\tau_1-\tau_2)d\tau d\tau_1d\tau_2\\=\int_{-\infty}^\infty \int_{-\infty}^\infty x(\tau_1)x^*(\tau_2)R_h(t+\tau_1-\tau_2) d\tau_1d\tau_2\\=\int_{-\infty}^\infty \int_{-\infty}^\infty x(\tau_1+\tau_2)x^*(\tau_2)R_h(t+\tau_1) d\tau_1d\tau_2\\=\int_{-\infty}^\infty R_h(t+\tau_1) \int_{-\infty}^\infty x(\tau_1+\tau_2)x^*(\tau_2) d\tau_2 d\tau_1\\=\int_{-\infty}^\infty R_h(t+\tau_1) \int_{-\infty}^\infty x(\tau_2)x^*(\tau_2-\tau_1) d\tau_2 d\tau_1\\=\int_{-\infty}^\infty R_h(t+\tau_1) R_x(-\tau_1) d\tau_1\\=R_x(t)*R_h(t)
}$$

  (2) $x(t)$ is stochastic and $h(t)$ is deterministic.

In this case, $y(t)$ is stochastic and we obtain:
$$R_y(t_1,t_2){=\Bbb E\{y(t_1)y^*(t_2)\}
\\=\Bbb E\left\{\int_{-\infty}^\infty\int_{-\infty}^\infty x(\tau_1)h(t_1-\tau_1)x^*(\tau_2)h^*(t_2-\tau_2)d\tau_1d\tau_2\right\}
\\=\int_{-\infty}^\infty\int_{-\infty}^\infty \Bbb E\{x(\tau_1)x^*(\tau_2)\}h(t_1-\tau_1)h^*(t_2-\tau_2)d\tau_1d\tau_2
\\=\int_{-\infty}^\infty\int_{-\infty}^\infty R_x(\tau_1,\tau_2)h(t_1-\tau_1)h^*(t_2-\tau_2)d\tau_1d\tau_2
\\=\int_{-\infty}^\infty\left[\int_{-\infty}^\infty R_x(\tau_1,\tau_2)h(t_1-\tau_1)d\tau_1\right]h^*(t_2-\tau_2)d\tau_2
\\=\int_{-\infty}^\infty\left[ R_x(t_1,\tau_2)*_{t_1}h(t_1)\right]h^*(t-\tau_2)d\tau_2
\\=\left[ R_x(t_1,t_2)*_{t_1}h(t_1)\right]*_{t_2}h^*(t_2)
}$$where $*_u$ means convolution w.r.t. $u$.

  (3) $h(t)$ is stochastic and $x(t)$ is deterministic.

This case is very similar to case (2). Since the system is LTI, we can reverse the impulse response with input signal to obtain $$R_y(t_1,t_2)=\left[ R_h(t_1,t_2)*_{t_1}x(t_1)\right]*_{t_2}x^*(t_2)$$

  (4) Both $x(t)$ and $h(t)$ are stochastic.

As the statistical sources of $x(t)$ and $h(t)$ are often independent, it is quite reasonable to assume that $x(t)$ and $h(t)$ are also independent. To use this fact, we start by:$$R_y(t_1,t_2){=\Bbb E\{y(t_1)y^*(t_2)\}
\\=\Bbb E\left\{\int_{-\infty}^\infty\int_{-\infty}^\infty x(\tau_1)h(t_1-\tau_1)x^*(\tau_2)h^*(t_2-\tau_2)d\tau_1d\tau_2\right\}
\\=\Bbb \int_{-\infty}^\infty\int_{-\infty}^\infty \Bbb E\left\{x(\tau_1)h(t_1-\tau_1)x^*(\tau_2)h^*(t_2-\tau_2)\right\}d\tau_1d\tau_2
\\=\Bbb \int_{-\infty}^\infty\int_{-\infty}^\infty \Bbb E\{x(\tau_1)x^*(\tau_2)\}\Bbb E\{ h(t_1-\tau_1)h^*(t_2-\tau_2)\}d\tau_1d\tau_2
\\=\Bbb \int_{-\infty}^\infty R_x(t_1,\tau_2)*_{t_1}R_h(t_1,t_2-\tau_2)d\tau_2
\\=R_x(t_1,t_2)*_{t_1}*_{t_2}R_h(t_1,t_2)
}$$
which means that a dual convolution must be performed: one on $t_1$ and one on $t_2$.


what is the correct answer for shifting and scaling this signal? The discrete signal is provided in the following table:

$$\begin{array}{|c|ccccccc|}\hline
 n& -3 & -2 & -1 & 0 & 1 & 2 & 3\\ \hline
x[n] & 3 &  2 &  1 & 0 & 1 & 2 & 3 \\ \hline
\end{array}$$

The book is asking to determine $x[3n-1]$. This involves shifting by one and then scaling by $3$. My answer to this problem is:

$$\begin{array}{|c|cc|}\hline
 3n-1 &    0 & 1\\ \hline
x[3n-1]&  1 & 2 \\ \hline
\end{array}$$

However the answer in the book is:

$$\begin{array}{|c|ccc|}\hline
 3n-1 &  -1 & 0 & 1\\ \hline
x[3n-1]& 0 & 1 & 2 \\ \hline
\end{array}$$

Which answer is correct?

The sequence $x[n]$ is only non-zero in the range $-3\le n\le 3$. So if you change the index to $3n-1$, in order to remain inside that region, $n$ can only take the values $n=0$ (giving $3n-1=-1$), and $n=1$ (giving $3n-1=2$). This gives

$$x[3n-1]=\begin{cases}x[-1]=1,&n=0\\x[2]=2,&n=1\\
0,&\text{otherwise}\end{cases}$$

So if the exercise is as you stated it, then your solution is correct, apart from the fact that it is not $3n-1$ but $n$ that takes the values $0$ and $1$.


Mathematics behind the creation of 3d face mask at runtime I'm new in Image processing and Machine learning area. I want to create 3d mask of face at runtime.
import numpy as np
import dlib

cap = cv2.VideoCapture(0)

detector = dlib.get_frontal_face_detector()
predictor = dlib.shape_predictor("shape_predictor_68_face_landmarks.dat")

while True:
    _, frame = cap.read()
    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    faces = detector(gray)
    for face in faces:
        x1 = face.left()
        y1 = face.top()
        x2 = face.right()
        y2 = face.bottom()

        landmarks = predictor(gray, face)

        for n in range(0, 68):
            x = landmarks.part(n).x
            y = landmarks.part(n).y
            cv2.circle(frame, (x, y), 3, (255, 0, 0), -1)

    cv2.imshow("Frame", frame)

    key = cv2.waitKey(1)
    if key == 27:
        break

I find mention code online, which find 68 point on face.
I want to know the Mathematics behind it. All I know that partial differentiation is used but I can't understand how. Also Some details about how to create 3d face mask at runtime will be appreciated.

Sorry for the naive answer, but why not consider reading the documentation of the library you have used?
http://dlib.net/
http://dlib.net/python/index.html#dlib.get_frontal_face_detector
http://dlib.net/python/index.html#dlib.shape_predictor
Quick googling also reveals numerous posts on how this library works and how it can be used
https://www.pyimagesearch.com/2019/12/16/training-a-custom-dlib-shape-predictor/
DLIB appears to be a machine learning / deep learning library, which means they would train a predictor like a deep neuronal network on some training data, and then use that network to place points on your images. The mathematics behind that is no different than any other neuronal network. Just find some course on machine learning / deep learning, and you will learn what you need.


Hardware implementation of Gaussian noise generator I am implementing a GNG using Verilog. The implementation is based on the IEEE paper "Hardware Architecture of a Gaussian Noise Generator Based on the Inversion Method". In this, i could not understand what exactly does the ROM_trans block needs to store. I understand it is a lookup table, but I do not know what to store inside it. Also how will I determine the coefficients that need to be stored in the ROM_coef block.

My friend, you dont need to break your head against a paper...

Just with any example of this list, the Blum Blum Shub Method:

$$x'_{{n+1}}=x_{n}^{'2}{\bmod M}$$

with $M=pq$, $p$, $q$ primes, with $x=\frac{1}{M}x' \approx U(0,1)$ gives you a fairly acceptable uniform $U(0,1)$ random number generator. You can use $p$=13331 and $q$=131 or whatever other.

And for converting the $U(0,1)$ Uniform Random Number into a $N(0,1)$ Gaussian Number Generator, as in here, just can use the Central Limit Theorem to generate N numbers each step, and make the sum (Bates Distribution). You can use $N$=15:

$$y'_{n}=\frac{1}{N}\sum_{i=1}^N x_{n}^i$$

with $y=\sqrt{12N}(y'-\frac{1}{2}) \approx N(0,1)$

For randomizing, just use a clock or cpu time seed as $x_0'$.


Where is the flaw in this derivation of the DTFT of the unit step sequence $u[n]$? This question is related to this other question of mine where I ask for derivations of the discrete-time Fourier transform (DTFT) of the unit step sequence $u[n]$. During my search for derivations I found one which is amazingly simple. I first saw it on page 138 of this book by B.A. Shenoi. I also came across it on mathematics.SE in this answer.

Since the argument is short and simple I will repeat it here for convenience.

  The unit step sequence can be written as
  $$u[n]=f[n]+\frac12\tag{1}$$
  with
  $$f[n]=\begin{cases}\frac12,\quad n\ge 0\\-\frac12,\quad n<0\end{cases}\tag{2}$$
  Obviously,
  $$f[n]-f[n-1]=\delta[n]\tag{3}$$
  Applying the DTFT on both sided of $(3)$ gives
  $$F(\omega)\left(1-e^{-j\omega}\right)=1\tag{4}$$
  where $F(\omega)$ is the DTFT of $f[n]$. From $(4)$ we get
  $$F(\omega)=\frac{1}{1-e^{-j\omega}}\tag{5}$$
  From $(5)$ and $(1)$ we get for the DTFT of $u[n]$
  $$U(\omega)=F(\omega)+\pi\delta(\omega)=\frac{1}{1-e^{-j\omega}}+\pi\delta(\omega),\quad -\pi\le\omega <\pi\tag{6}$$
  where I've used $\text{DTFT}\{1\}=2\pi\delta(\omega)$, $-\pi\le\omega <\pi$.

Eq. $(6)$ for the DTFT of $u[n]$ is no doubt correct. However, the derivation is flawed.

The question is: find and explain the flaw in above derivation.

Please prepend your answer with the spoiler tag >!.

There are infinitely many signals that make the following equality hold: $$y[n]-y[n-1]=\delta[n] \qquad (1)$$
The only thing that matters is that $y[0]-y[-1]=1$, and then the rest of the coefficients of $y$ can be determined under the restriction that Eq. $(1)$ states (i.e. the substraction of consecutive samples must be $0$ for $n\neq 0$).
In other words, Eq. $(1)$ will be achieved by any signal $y[n]$ such that $$y[0]=y[-1]+1 \land y[n]=y[n-1] \ \forall n\neq0$$
Another way to see this is that any function that is basically $u[n]$ with an offset (a constant value added) will satisfy $(1)$. This explains the statement made by robert bristow-johnson in his answer: differentiators destroy this information (such as taking a derivative in continuous time destroys evidence of any constant value in the original function).

To sum up, I believe that the proof is flawed because the procedure followed could use any function of the form $u[n]+C$ with $C\in\mathbb{R}$, and this would lead to many functions having the same Fourier transform, which is indeed wrong as the Fourier transform is a bijection. Maybe the author deliberately decided to ignore anything related to DC values, conscious that in order to show that $F(\omega)$ is the DTFT of $f[n]$ he would need the accumulation property (whose most popular proof is derived from the DTFT of the unit step - ergo, a pretty circular proof). The proof is not strictly wrong, as everything it states (the formulae for $F(\omega)$ and $U(\omega)$, the decomposition of the unit step, the difference equation) is true, but it would require the accumulation property to show why $F(\omega)$ doesn't have any Dirac deltas.


Odd sampling for symmetric functions? I have just started my first proper course in DSP, and am a little confused about a statement in my textbook.  The book states:

"Notice that when the Nyqust-frequency is defined as $\pi/\Delta t$, we will with $\omega_m = 2 \pi m/N \Delta t$ where $m = 0, \pm 1, \pm2, ...$ have $\omega_m \neq \pi/ \Delta t$, unless $N$ is even.  If we are working with symmetric functions, it can be wise to use an odd number of sampling points."

I see the logic in the first sentence here, but don't understand how it follows from this that it is wise to use an odd number of sampling points for a symmetric function. This statement is not elaborated upon further. Can anyone explain this logic to me?  Why it is wise to use an odd number of sampling points for a symmetric function?  If anyone could explain this, I would greatly appreciate it!

It is because it is symmetric about m = 0. So if you have an odd number of points then you can also sample at +/-1, +/-2, and so on. Thus each of the points has it's counterpart except for zero. This only works with an odd number of points.


DWT architecture using filter bank I am studying this paper.

In one of the figures in the paper (Fig 5.b and Fig 5.c) DWT architecture is given using db2 and filter bank.

I don't understand how Lo_D and Hi_D have size half of input image size? If we see the figures there is no downsampler. Then how the size of the Lo_D and Hi_D have half the size of image (either row wise or column wise)?

Please explain me how single level DWT is implemented using filter banks? How the size of the decomposition values is half of the input (i.e. image) when filter banks are used to implement DWT?

I understand they  implement the convolution by Lo and Hi db2 filters given after the sentence:

  Rather, we propose to use db2 filters. The filter coefficients for the
  db2 filter are represented as

I have not seen in the paper a claim that Lo_D and Hi_D in Figures 5b and 5c should be half-size. They are low-pass and high-pass branches. However, they should be followed by a downsampling to implement a  one level-DWT, or at least something similar.

There are different ways for  implementing a one level of a standard critically sampled DWT: direct, lattice, lifting, polyphase, Fourier. They should end up all with the  same number of samples as the original signal (away from the borders).


What's the sampling frequency when synthetically generating data? Say I am using matlab's random number generator to get an array of random numbers, which I will then, for testing purposes, find the spectral density etc. For that purpose I need to know the sampling freq., but these have been generated synthetically so I am not sure how the sampling theorem applies here.

The sampling rate is arbitrary- you choose it.

Usually people are simulating a real situation, like, for example, audio from a CD.  Then you would choose the sample rate to be 44.1 kHz.


IFFT Concept in OFDM can any one explain how actually IFFT generating sub carriers in OFDM.... when we are simulating OFDM transmitter in matlab, we are taking input data stream, encoding it and modulating the encoded data. then we are converting this modulated serial data stream to parallel data and giving it to the IFFT block. my problem is how sub carriers are generating in IFFT??

The first thing you must change for understanding why IFFT is used for OFDM is to forget about time and frequency domains.

As the name states OFDM is a technique about frequency domain multiplexing. When you think in FDM you'll have several channels modulated with a guard interval between then to prevent interference. If the carriers are orthogonal you'll not have interference, but without the waste of band guard.

But how you can have this?

One way is to use several modulators, as in FDM, but with the correct frequencies. This is hard to accomplish.

If you will implement digitally the filter bank,and this is the correct way to view IFFT on OFDM systems: a transmultiplexer(TMUX), is equivalent to an IDFT. This is why you can use IFFT as a fast implementation for OFDM. You can read one of the books about multi-carrier communication for a better explanation. This paper
http://wsl.stanford.edu/~ee359/ofdm_bingham.pdf is also one that I like about multi-carrier modulation.


What is the time-integration property in the Fourier series analysis? In the continuous Fourier series properties for a periodic continuous-time signal, we have time-integration property:

$$
\int_{-\infty}^t x(\alpha)d\alpha \leftrightarrow \frac{a_k}{jk\omega_0}
$$

where $a_k$ is the Fourier series coefficients of the signal $x(t)$. Now I am confused. Is the above integral equivalent to $\int x(t)dt$ ?

For instance, the Fourier series coefficients of $\cos(\omega_0 t)$ are: $a_1=a_{-1}=\frac{1}{2}$, other $a_{k}=0$.

and the Fourier series coefficients of  $\sin(\omega_0 t)$ are: $b_1=\frac{1}{2j}, b_{-1}=\frac{-1}{2j}$, other $b_{k}=0$.

and we have: $\int \cos(\omega_0t)dt=\frac{1}{\omega_0}\sin(\omega_0t)$

Then:

$$
\textrm{F.S.}\left\{\int \cos(\omega_0t)dt\right\} = \frac{a_k}{jk\omega_0} = \textrm{F.S.}\left\{\frac{1}{\omega_0}\sin(\omega_0t)\right\}= \frac{1}{\omega_0}b_k
$$

Thus:

$$
b_k=\frac{a_k}{jk}
$$

as we expected.

Let's consider the complex form of the Fourier series of a $T$-periodic function $x(t)$:

$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{jk\omega_0t}\tag{1}$$

with $\omega_0=2\pi/T$. The integration property says that the Fourier series of

$$y(t)=\int x(t)dt\tag{2}$$

is given by

$$y(t)=\sum_{k=-\infty}^{\infty}b_ke^{jk\omega_0t}\tag{3}$$

with

$$b_k=\frac{a_k}{jk\omega_0},\quad k\neq 0\tag{4}$$

where we require that $a_0=0$.

Note that in $(2)$ I used the indefinite integral, and I chose the integration constant to be zero, such that $b_0=0$. In some texts you may find the definite integral instead of the indefinite integral:

$$y(t)=\int_{-\infty}^{t} x(t)dt\tag{5}$$

The problem with $(5)$ is that for certain functions $x(t)$ the definite integral doesn't exist with $-\infty$ as the lower bound. Note that the lower bound in $(5)$ corresponds to the integration constant of the indefinite integral $(2)$. This becomes obvious if we consider the relation between the indefinite and the definite integral:

$$\int x(t)dt=X(t)+C\tag{6}$$

$$\int_{t_0}^tx(\tau)d\tau=X(t)-X(t_0)\tag{7}$$


Why should I zero-pad a signal before taking the Fourier transform? In an answer to a previous question, it was stated that one should

  zero-pad the input signals (add zeros to the end so that at least half of the wave is "blank")

What's the reason for this?

Zero padding allows one to use a longer FFT, which will produce a longer FFT result vector.

A longer FFT result has more frequency bins that are more closely spaced in frequency.  But they will be essentially providing the same result as a high quality Sinc interpolation of a shorter non-zero-padded FFT of the original data.

This might result in a smoother looking spectrum when plotted without further interpolation.

Although this interpolation won't help with resolving or the resolution of and/or between adjacent or nearby frequencies, it might make it easier to visually resolve the peak of a single isolated frequency that does not have any significant adjacent signals or noise in the spectrum.  Statistically, the higher density of FFT result bins will probably make it more likely that the peak magnitude bin is closer to the frequency of a random isolated input frequency sinusoid, and without further interpolation (parabolic, et.al.).

But, essentially, zero padding before a DFT/FFT is a computationally efficient method of interpolating a large number of points.

Zero-padding for cross-correlation, auto-correlation, or convolution filtering  is used to not mix convolution results (due to circular convolution).  The full result of a linear convolution is longer than either of the two input vectors.  If you don't provide a place to put the end of this longer convolution result, FFT fast convolution will just mix it in with and cruft up your desired result.  Zero-padding provides a bunch zeros into which to mix the longer result.  And it's far far easier to un-mix something that has only been mixed/summed with a vector of zeros.


Is it wise to use mathematical morphology and Gabor filter to extract blood vessels in the same algorithm? I'm just curious, is it really wise to use mathematical morphology and Gabor filtering in the same code to extract blood vessels?

Like in this case, I use morphology to get the vessels and using Gabor filter to further enhancing brightness of my vessel.

In a word, yes.

Mathematical morphology and harmonic analysis (Gabor for instance) are two families of tools, relying on different sets of assumptions, pertaining to image formation models. Harmonic analysis dwells in vector spaces, addressing better linear or superposition  effects. Mathematical morphology, built on lattices, could be more adapted to connected objects, discretized levels, non-linear effects. Both types of contents coexist in images, so mathematical morphology and  harmonic analysis can be complementary tools for complicated tasks.

Some tools reside at the intersection, like morphological wavelets, some links relate linear filters and structuring elements (logarithmically, through the Cramer transform, see fro instance An Explanation for the Logarithmic Connection between Linear and Morphological System Theory). However, it remains wise to use the right set of tools when it is at its best.

For instance, in the most simple signal processing applications, where information is degraded by both  impulse and Gaussian  noise, it is common to first use a non-linear median filter, and a linear low-pass one.

Finally, let me mention a recent thesis related to your topic: RORPO : A morphological framework for curvilinear structure analysis; Application to the filtering and segmentation of blood vessels, where mathematical morphology and harmonic analysis are discussed.


Denoise EEG signal by using Daubechies function I have an EEG signal and it contains eye blink artifacts. I read some references and know that it is possible to detect eye blinks and remove them by using wavelet transforms, but I don't know that how do it.

How do I detect eye blinks? After transforming the EEG signal into wavelet coefficients, what should I do and which level of Daubechies can be used to do that?

You might simply use soft/hard thresholding. This operation is quite standard and called Wavelet Denoising.
Here are some resources:

http://www.stthomas.edu/mathematics/pdfs/MSAD/Denoising%20via%20Wavele.pdf
http://scholar.lib.vt.edu/theses/available/etd-12062002-152858/unrestricted/Chapter4.pdf
http://www.ee.iitb.ac.in/~icvgip/PAPERS/202.pdf

An overview of noise removal for EEG:
http://www.ijarcce.com/upload/february/3_A%20survey%20on%20different.pdf

Some source code on image denoising via VisuShrink:
http://www.mathworks.com/matlabcentral/fileexchange/43996-image-denoising-using-visushrink


Refining lines in an image with openCV In my project I am extracting the lines of an image with openCV2 by firstly creating a skeleton: I convert it to a binary image, invert it (to make it work with morphological operations) and then perform a bunch of morphological operations. Then I apply the canny edge detector on the skeleton and perform a Hough transformation.
Now the important part: In order to remove the gaps I perform a closing operation. My problem now is that there are still some gaps and fragments and I can't apply the closing operation one more time because its idempotent (only applicable once).
I want to refine the lines by closing all gaps and possibly thicken the lines. How may I achieve this?
Here some images to illustrate my problem.
Original image

Processed image

See those tiny lines which protrude from the edges? Or the interruptions in a straight line?
For the gaps in the lines I thought about looping through the matrix and if the points left and right of a pixel are 1 (drawn), then that particular pixel should also be set to 1. Would that be a good solution?

If you are trying to remove the gaps during hough transform, you can set minLineLength(Minimum length of line. Line segments shorter than this are rejected) and maxLineGap (Maximum allowed gap between line segments to treat them as single line) in cv2.HoughLinesP()

EDIT

Some results from the morphological operations:


Difference between bandpass filter and using fast fourier transform to extract power of frequency bands In signal processing, I have heard of two terms:

Using band-pass filter to extract some frequency bands.
Using Discrete/Fast Fourier Transform (DFT/FFT) to extract some frequency bands.

Are these two terms similar or different?

How and when we use each of these methods?

Thanks in advance

Fast Fourier transform does not extract any frequency bands. It only shows the frequency content of a given signal. But while applying FFT, one should be careful about choosing the sampling frequency. As an example, if a signal contains a frequency range of $0-100\,\text{Hz}$ and $f_{max}=100\,\text{Hz}$ $$f_{max} = \dfrac{F_s}{2}$$ and the sampling frequency $F_s \geq 2\cdot f_{max}$ of the signal should be greater than $200\,\text{Hz}$ (according to Nyquist frequency rule). Then FFT will show all the frequencies present in the signal i.e. $0-100\,\text{Hz}$.

If the sampling frequency is less than $<200\,\text{Hz}$ then aliasing will occur and will provide an aliased frequency representation.

On the other hand, band-pass filters eliminate some frequency content of the signal. As an example say one signal has $0-100\,\text{Hz}$ frequency components but you are interested in exploring the frequencies of range $10-30\,\text{Hz}$, then one should use a bandpass filter.

As a system identification researcher, I can provide an example where FFT and bandpass filters are used. Say a signal has its main frequencies in the $10-30\,\text{Hz}$ range. But the signal is acquired with a sampling frequency of $200\,\text{Hz}$.

Performing FFT will show frequency components from $0-100\,\text{Hz}$. Say there are measurement noises all over the frequency ranges. I use a bandpass filter to get rid of that measurement noise.


Question about eigendecomposition, signal subspace and their properties Lately, I am readig a paper titled
A Subspace Method for Estimating Sensor Gains and Phases.

In it, it is mentioned:

There are $m$ sensors in the array, $n$ known  narrowband far-field signals($m \gt n$, equal power). $R$ is covariance matrix of received signal, $\Gamma \in \mathbb{C}^{m \times m}$ is the diagonal matrix of unknown  sensor gains and phases, which this paper tried to solve. $A\in \mathbb{C}^{m \times n}$ is the steering matrix. $E_s$ are eigenvectors associated with principal eigenvalues. $^\dagger$ represents conjugate transpose.

It's easy to know $E_s \in \mathbb{C}^{m \times n}$ and $rank(E_s)=n.$

My first question: I can't figure out why the eq(6) holds. Did I miss some properties of the signal subspace $E_sE_s^H?$

Then the equation was transformed into eq(9),

where $D_i$ is each column of $A$(in diagonal form) and $v$ is column vector formed form diagonal elements in $\Gamma$.

My second question: Why $W_i=D_i^HE_sE_s^HD_i$ has n unity eigenvalues and m-n zero eigenvalues?(underlined in blue).

Thank for your patience to read my questions!

My first swing at the answer had some very incorrect claims.

I do not have access to the article, so I am inferring some things from the portion posted in the question.

NOTA BENE: My arguments assume that the eigenvectors of $\mathbf{R}$ are arranged so that the first $n$ belong to the signal subspace and that the last $m-n$ belong to the noise subspace. That is how the formulas appear so clean.

Since $\mathbf{R}$ is a covariance matrix, it is a positive-definite Hermitian matrix. This ensures that

all of its eigenvalues are positive,
it is unitarily diagonalizable.

That last property means that there is a unitary matrix $\mathbf{U}$, whose columns are orthonormal eigenvectors of $\mathbf{R}$, such that $\mathbf{U}^{\dagger}\mathbf{R}\mathbf{U}$ is equal to a diagonal matrix that has $\mathbf{R}$'s eigenvalues on its diagonal:
\begin{eqnarray}
\mathbf{U}^{\dagger}\mathbf{R}\mathbf{U} = \textrm{diag}(\lambda_1,\ldots,\lambda_n,\underbrace{\sigma^2,\ldots,\sigma^2}_{\textrm{$m-n$ $\sigma^2$s}}).
\end{eqnarray}

Because the article mentions $\mathbf{E}_s^{\dagger}$ and $\mathbf{E}_n^{\dagger}$ and not $\mathbf{E}_s^{-1}$ and $\mathbf{E}_n^{-1}$ in the eigenvalue-eigenvector decomposition of $\mathbf{R}$, I suspect that the unitary diagonzaling matrix of $\mathbf{R}$ is $\mathbf{E}_s + \mathbf{E}_n$, where $\mathbf{E}_s^{\dagger}\mathbf{E}_n = \mathbf{O}$, an all-zero matrix. I will demonstrate using the decomposition given in the article.
\begin{eqnarray}
&& (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{\mathbf{R}}(\mathbf{E}_s + \mathbf{E}_n)\\
&=& (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{(\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n\mathbf{E}_n^{\dagger})}(\mathbf{E}_s + \mathbf{E}_n)\\
&=& (\mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger})(\mathbf{E}_s + \mathbf{E}_n),
\end{eqnarray}
where I have already dropped $\mathbf{E}_s^{\dagger}\mathbf{E}_n$ and $\mathbf{E}_n^{\dagger}\mathbf{E}_s$, both of which are all-zero matrices. The next step is
\begin{equation}
\mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger}\mathbf{E}_s + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger}\mathbf{E}_n.
\end{equation}

The OP states that $\mathbf{E}_s$ is $m\times n$ and has rank $n$, where $n < m$. But I suspect (hope) that $\mathbf{E}_s$ is $m\times m$ with rank $n$ and that and $\mathbf{E}_n$ is $m\times m$ with rank $m-n$, with the following structures:
\begin{eqnarray}
\mathbf{E}_s &=& (\mathbf{u}_1\cdots\mathbf{u}_n\underbrace{\mathbf{0}\cdots\mathbf{0}}_{m-n}),\\
\mathbf{E}_n &=& (\underbrace{\mathbf{0}\cdots\mathbf{0}}_{n}\mathbf{u}_{n+1}\cdots\mathbf{u}_m)
\end{eqnarray}
That is, $\mathbf{E}_s$ should have $n$ orthonormal columns followed by $m-n$ all-zero columns, while $\mathbf{E}_n$ starts with $n$ all-zero columns and ends with $m-n$ orthonormal columns.

If these conditions were true, then we would have
\begin{eqnarray}
\mathbf{E}_s^{\dagger}\mathbf{E}_s &=& \textrm{diag}(\underbrace{1,\ldots,1}_{n},\underbrace{0,\ldots,0}_{m-n}),\\
\mathbf{E}_n^{\dagger}\mathbf{E}_n &=& \textrm{diag}(\underbrace{0,\ldots,0}_{n},\underbrace{1,\ldots,1}_{m-n})
\end{eqnarray}
and
\begin{eqnarray}
&& (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{\mathbf{R}}(\mathbf{E}_s + \mathbf{E}_n)\\
&=& \mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger}\mathbf{E}_s + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger}\mathbf{E}_n\\
&=& \textrm{diag}(\lambda_1,\ldots,\lambda_n,\underbrace{\sigma^2,\ldots,\sigma^2}_{m-n})
\end{eqnarray}

Even if $\mathbf{E}_s$ is $m\times n$ and $\mathbf{E}_n$ is $m\times(m-n)$, we still have many of the same properties as long as
\begin{eqnarray}
\mathbf{E}_s &=& (\mathbf{u}_1\cdots\mathbf{u}_n),\\
\mathbf{E}_n &=& (\mathbf{u}_{n+1}\cdots\mathbf{u}_m).
\end{eqnarray}
In particular, we have
\begin{equation}
\mathbf{E}_s^{\dagger}\mathbf{E}_n = \mathbf{O}_{n\times(m-n)}
\end{equation}
and
\begin{equation}
\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n\mathbf{E}_n^{\dagger}
=
\sum_{k=1}^{n}\lambda_k\mathbf{u}_k\mathbf{u}_k^{\dagger}
+
\sum_{k=n+1}^{m}\sigma^2\mathbf{u}_k\mathbf{u}_k^{\dagger},
\end{equation}
which is a spectral decomposition of a Hermitian matrix such as $\mathbf{R}$.


Arbitrary-Response Filter Design I want to design a FIR filter by using  'cfirpm' function , which can i use the group delay and grid frequency,
my code is :

N=30 ;
f=[linspace(-1,-.5,120),linspace(-0.4,0.7,60),linspace(0.8,1,90)];
d=[zeros(1,120),ones(1,60),zeros(1,90)].*exp(-j*pi*f*N-1/2);
b=cfirpm(N,d,@lowpass);
fvtool(b,1,'OverlayedAnalysis','phase');

But i got the following error in matlab,

Error using cfirpm
Band edges must be monotonically increasing.

How can i solve this problem? or,how can i add the group delay in this function?

Thank you in advance.

For a low-pass filter it is sufficient to specify the band edges, as shown in this example on the mathworks website. From your code I suppose that your band edges are

f=[-1,-0.5,-0.4,0.7,0.8,1]

Note that you wrote -0.7 and -0.8 (i.e., with a negative sign), which I assume is a mistake.

Furthermore, think about the desired phase or group delay response. The phase factor exp(-j*pi*f*N-1/2) does not make much sense, an it is probably not what you mean. If you want a linear phase (with delay $N/2$, where $N$ is the filter order) there is no need to use cfirpm, and you could better use firpm.


EEG, PSI interpretation My colleague and I started exploring EEG processing and, especially, functional connectivity. We decided to try the Phase Slope Index (PSI) and conducted a simple experiment as follows: it consists of three steps, research participant should raise his left hand (1), then his right one (2), then both (3). The description is simplified a bit but shows what the experiment was. We expected to see the information flow between the pre-motor cortex and primary motor cortex but we seemed to fail. The data was pre-processed (filtered at [0.25; 40] Hz, re-referenced to the average electrode, segmented into three pieces according to three stages of experiment; also blink artefacts were removed, so were motion ones). Then we computed the all-to-all PSI matrix and what we have had is depicted on the pictures (if needed and more appropriate, I could add the graph representation).   All in all, we haven't seemed to reach what we aimed and now we're trying to figure out whether the fail refers to the interpretation, or experiment design, or both, or something else. My questions are:

Can one guess what stage of experiment every matrix refers to? If so, how? I mean, do they at least look like the result of the experiment described? If they do, how can one interpret them?
What could be an experiment that would show the noticeable PSI for two areas given beforehand? I mean, we had two areas given, pre-motor cortex and primary motor cortex, and we supposed that we would find a synchronization between them in such a motion experiment; something's gone wrong, though. If the problem is in the design of experiment, what design would be better?

Can one guess what stage of experiment every matrix refers to? If so, how? I mean, do they at least look like the result of the experiment described? If they do, how can one interpret them?

I do not think that I can. But here is a guess that might give you some insight: The most "populated" matrices are those for a single hand up in the air (whether left or right).

Why?

Functional connectivity metrics (of which there are so many by now) assess some of "dependence" or "correlation". There is a higher chance for any two pairs of waveforms to have more correlated signals when the diversity of the generated signals is low. With two hands being raised, the brain generates a large set of impulses to drive all the necessary muscles. Consequently, the chance of them having a high correlation is low.

If you notice, a key characteristic of functional connectivity estimation in related papers is some form of averaging, whether temporal, spatial or across and within subjects.

Whether it is Magnetoencephalography (higher spatial resolution than Electroencephalography), or of course, fMRI or EEG, there is an intermediate stage that "summarises" observations.

In fMRI you sum voxel traces around anatomical areas indicated by an atlas, in MEG you sum contributions of elements around broad areas (simplest method) OR you do some form of source localisation and then look at the connectivity of the localised sources. Similar workflow applies in EEG too (e.g. LORETTA as above).

This is a required step to increase the Signal to Noise Ratio (SNR) and look at (hopefully) robust correlations. Along this idea are also methods that examine functional connectivity in specific bands of brain activity

Therefore, I doubt that you will get anything meaningful by looking at each individual channel without some form of spatial or temporal filtering / averaging.

This however, is half of the story. The other thing you can do is, of course, to change the task.

  What could be an experiment that would show the noticeable PSI for two areas given beforehand? I mean, we had two areas given, pre-motor cortex and primary motor cortex, and we supposed that we would find a synchronization between them in such a motion experiment; something's gone wrong, though. If the problem is in the design of experiment, what design would be better?

You will definitely see differences in the default mode network between two radically different tasks. For example, at rest and at "work".

To reduce your post-processing workload, try to use tasks that are not likely to generate a lot of artifacts. For example: Resting state and "solve a mathematical problem".

This design will give you something to compare "activity" against and help in reducing spurious correlations further.

Finally, be aware of filtering. The more your filter a waveform, the more it "looks like the filter" and therefore, the more chances there are to have more increased correlations.

Hope this helps.


why there is no 0 on the right side? (matlab, low pass filter) I took even number of sample (in this case 6) from sinc, hamming window, and multiple of those two.
then after fourier transform, frequency-domain variables have zero only left side, in this case its x is -pi.
I'm wondering why there are no zeros in right side (where x = pi), and has positive value.

here's my work using fourier transfrom

Answer: Your derivation is absolutely correct and You will see 0 for both $\omega = -\pi$ and $\omega = \pi$, if you use fvtool() function of MATLAB instead of fft(). The reason is simple : FFT does not calculate values of $H(e^{j\omega})$ at continuous $\omega \in [-\pi, \pi]$, but it calculates DFT and for that the digital frequency resolution is $\omega = \frac{2\pi k}{N}$. So, for $N=6$, it will compute values of $H(e^{j\omega})$ at $\omega = 0, \frac{2\pi}{6}, \frac{4\pi}{6}, \frac{6\pi}{6}, \frac{8\pi}{6} \ and \ \frac{10\pi}{6}$. See that $\omega = \pi$ is occurring only at $k=3$, and hence you are seeing $0$ only once.

For checking this, you can take any symmetrical h, like $h = {1,2,4,4,2,1}$.
Then use the following MATLAB code to plot its DTFT for $\omega \in [-\pi, \pi]$

  fvtool(h,1);

This will plot a Single sided DTFT from $\omega \in [0, \pi]$. You can then change from "Analysis" option to plot for frequency range $[-\pi, \pi]$.

You will get following plot for the example I have taken:

Check that the values are $0$ at both $\omega = -\pi$ and $\omega = \pi$.

And this will happen for any symmetric values, not any particular window function.


Fusing or combining signals from multiple sensors I've become very interested in the problem of fusing multiple low-cost sensors' outputs and trying to combine these outputs in such a way as to rival or exceed the output of a single high-quality / expensive sensor.

Applications could range from cubesats to antenna or microphone areays.

Specifically, it seems intuitively that one should be able to do better than simple "averaging" to, say, increase SNR for detecting faint signals or perhaps even increasing effective receiver bandwidth (i.e. combining multiple 3 MHz receiver outputs and achieving significantly greater BW without increasing sample rate).

I'm really not that smart on DSP, though, and would very much appreciate any pointers, theoretical limitations, suggested algorithms, resources, or online courses that could help me out.

Sensors fusion is a quite wide subtopic in signal processing. Apart from sensor fusion, it appears under different guises named for instance  diversity enhancement, multisensor combination, ensemble averaging, to name a few.

If the sensors' outputs $y_i$ are in phase or synchronized (a sensitive hypothesis since your are talking about averaging), but possibly with noise $n_i$ or different amplitude levels (the $k_i$ factor):

$$y_i = k_i x+n_i$$

there is a quite old reference in Signal extraction from multiple noisy sensors, S. Palit, 1997. I haven't yet implemented it, and it is not quite cited, but I still have the feeling it deserves a closer look. In the same genre, a lot of recent works follow the track of "ensemble averaging", and the notion of "low-rank" decomposition sounds appealing: put your outputs in a matrix $X$. If there were no noise, the columns would be multiple of the single signal $x$. Hence, is $x$ is not zero, the matrix $X$ is rank one. The challenge, with noise, is to find a low-rank approximation of $X$, separating the signal $x$, noises $n_i$ and potentially localized disturbances, outliers with a sparse footprint. One of these methods is known as GoDec: Randomized Low-rank & Sparse Matrix Decomposition in Noisy Case:

  Low-rank and sparse structures have been profoundly studied in matrix
  completion and compressed sensing. In this paper, we develop ``Go
  Decomposition'' (GoDec) to efficiently and robustly estimate the
  low-rank part $L$ and the sparse part $S$ of a matrix $X=L+S+G$ with
  noise $G$. GoDec alternatively assigns the low-rank approximation of
  $X-S$ to $L$ and the sparse approximation of $X-L$ to $S$. The
  algorithm can be significantly accelerated by bilateral random
  projections (BRP). We also propose GoDec for matrix completion as an
  important variant. We prove that the objective value $\|X-L-S\|_F^2$
  converges to a local minimum, while $L$ and $S$ linearly converge to
  local optimums. Theoretically, we analyze the influence of $L$, $S$
  and $G$ to the asymptotic/convergence speeds in order to discover the
  robustness of GoDec. Empirical studies suggest the efficiency,
  robustness and effectiveness of GoDec comparing with representative
  matrix decomposition and completion tools, e.g., Robust PCA and
  OptSpace. GoDec can be extended to solve multi-label learning problem
  by decomposing the multi-label data into the sum of several low-rank
  part and a sparse residual, where each low-rank part corresponds to
  the mapping of a particular label in the feature space. Then
  prediction can be obtained by finding the group sparse representations
  of a new instance on the subspaces defined by the low-rank parts.

If the outputs are not time-aligned, or with jitter, then you can explore convolutive blind source separation, superresolution (for BW extension), etc. but this margin is too thin to draw a complete panorama.


Measure the SNR of a signal in the frequency domain? How do you measure the SNR of a sinewave in the frequency domain? (Assuming no filtering.)

Suppose: $$x(n) = s(n) + n(n)$$

Create the signal (nCycles of a sinewave) in MATLAB...:

sampleRate = 1024;
f0 = sampleRate/8;

nCycles = 11;
time = 0:1/sampleRate:nCycles/f0;

signal = sin(2*pi*f0*time);

Create the noise and scale it by the desired SNR...

SNR = 10;
noise = randn(size(signal));
% scale the noise to obtain the desired SNR
noise = noise / norm(noise) * norm(signal) / 10^(SNR/20);

x = signal + noise;

Calculate the resulting SNR (in the time domain):

actualSNR = 20*log10(norm(signal)/norm(x - signal));
disp(['SNR = ',num2str(actualSNR),'  dB'])

Plot the noisy signal in time and frequency

subplot(2,1,1)
% plot the signal in time
plot(time,signal)
hold on; grid on
plot(time,x,'r')

% plot the noisy signal in frequency
NFFT = 4096;
X = fftshift(fft(x,NFFT));
X = X/max(abs(X));
f = sampleRate/2*linspace(-1,1,NFFT);
subplot(2,1,2)

plot(f,20*log10(abs(X)))
grid on;

It can be a little tricky.  If your sine wave happens to fall at an FFT bin center, things are a little easier.  If your sine wave happens to not fall at a bin center, you have to consider spectral leakage.  You'll also need to consider how your window function affects your signal.

A simple technique to estimate the signal power would be to sum the squared magnitudes of the FFT bins with most of the signal's energy (i.e. the largest bins).  To get the noise energy, sum the squared magnitude of all of the other bins.  This will get you close.


Calculating frequency and damping ratio from transfer function given eigenvalues I have the following standard transfer function for a damped linear oscillator:

$$G(s) = \dfrac{\omega_0^2}{s^2 + 2\zeta\omega_0s + \omega_0^2}$$

Now I have two eigen values at locations $-100 \pm 100i$.
I want to calculate the damping coefficient, $\zeta$, and natural frequency $\omega_0$.

I thought that the characteristic polynomial of matrix $A$ from the linear system in state space form is the denominator of $G(s)$. However, when solving $s^2 + 2\zeta\omega_0s + \omega_0^2=0$ for $s = -100 - 100i$, I get an equation with two unknowns right? What theoretical point am I missing here? how should I go about this?

Thanks in advance

If I understand you correctly your system has two complex conjugate poles at $p=-100+100i$ and $p^*=-100-100i$ (where $^*$ denotes complex conjugation). Consequently, your denominator polynomial can be written as
$$\begin{align}(s-p)(s-p^*)&=s^2-(p+p^*)s+|p|^2\\&=s^2-2\text{Re}\{p\}s+|p|^2\\&=s^2-2|p|\cos(\phi)s+|p|^2\end{align}\tag{1}$$
where $\phi$ is the pole angle: $p=|p|e^{i\phi}$. For the given poles you have $\phi=3\pi/4$. Comparing $(1)$ with the denominator in your question immediately gives $\omega_0=|p|$ and $\zeta=-\cos\phi$.


Effect of nonlinear transformation on digital signal Currently, I know that by passing a band-limited baseband digital signal through a nonlinear system results in the expansion of the original signal's bandwidth and creates harmonics.

Therefore, I think that I should upsample the original signal to a certain degree in order to avoid any higher frequency components from folding down to the lower frequency band after the conversion.

My question:

How should I determine the size of the expanded bandwidth? Does it have something to do with the highest harmonic visible above the noise floor?
Are there any ways to suppress the harmonics?
Do I need to pass the nonlinearly transformed signal through some kind of filter or upsample it further if I wish to feed it to a DAC?

Your intuition is correct. When simulating a discrete system, the sampling frequency needs to be high enough that you meet the Nyquist sampling criterion at all times, or you'll have folding.

A method I've used to determine the sampling frequency is to assume the signals are continuous, and find the Fourier transform of the system output. Usually you have an infinite number of harmonics, but you can usually identify a threshold after which you don't care about the harmonics folding. What the threshold is depends on your application and how much folding can you tolerate. Usually, harmonics 30 dB below your signal should pose no problems.
You can't suppress the harmonics at the nonlinear transformer's output, but you can filter them afterwards. Note that this will not help with the sampling frequency, since the harmonics are present between the transformer and the filter.
You need to make sure the signal is sampled at a frequency the DAC accepts. You can achieve that by resampling the signal right before the DAC.


Can the Walsh Hadamard transform be computed for non power of 2 sizes? Can the walsh hadamard transform be calculated for odd image block sizes such as 5x5 or 7x7? Most of the examples I've seen are for 4x4 and 8x8?

I fear it probably can't from the description I read on Wikipedia ( though I'm still trying to fully digest that page).

The Walsh-Hadamard transform requires a Hadamard matrix -- which has entries $\pm 1$ and
whose rows are orthogonal vectors.  A $n \times n$ Hadamard matrix is known to exist
for $n = 2$. Larger Hadamard matrices can exist only if $n$ is a multiple of $4$, though
it is not known if there is a Hadamard matrix for every multiple of $4$. However,
there is a recursive construction that can be used to construct a $2^m \times 2^m$
Hadamard matrix from a $2^{m-1}\times 2^{m-1}$ Hadamard matrix, and this structure
allows for the use of the Fast Hadamard Transform algorithm, which reduces the
computational cost from $2^{2m}$ additions and subtractions to $m2^m$ additions and
subtractions (just like the $N$-FFT reduces the number of multiplications from
$O(N^2)$ to $O(N\log N)$.

With this as background, the answer is Yes, you can use a Walsh-Hadamard transform
of length $N$ if you can find a $N\times N$ Hadamard matrix, but your choices
for $N$ are necessarily restricted. Also, fast algorithms may not exist
for your choice of $N$, though some speed-up is usually possible.

Note that Walsh-Hadamard transforms (WHTs) do not support cyclic convolutions, but
they do support what is sometimes called Poisson convolution. If $H$ denotes
a $2^m\times 2^m$
Hadamard matrix and $\mathbf x$ and $\mathbf y$ are vectors of length $2^m$ with \WHTs $\mathbf xH$ and $\mathbf yH$,
then the inverse WHT of the term-by-term multiplication of the entries
in $\mathbf xH$ and $\mathbf yH$ can be described as follows.  Suppose that the entries in
$\mathbf x$ and $\mathbf y$ etc are indexed not by integers $0$ through $2^m-1$ but rather
by the $m$-bit representations of these numbers.  Thus, we talk not of
$x[k]$ as the $k$-th entry in $\mathbf x$ but rather of $x[\mathbf k]$ where
$\mathbf k$ is the $m$-bit representation of $k$.  Then, the iWHT of the term-by-term
product of the entries in $\mathbf xH$ and $\mathbf yH$ has $\mathbf k$-th entry
$$\sum_{\mathbf i} x[\mathbf i]y[\mathbf k\oplus \mathbf i]$$
which is eerily reminiscent of
$$\sum_i x[i]y[k-i]$$
and even more so if one notes that modulo two, addition and subtraction
are the same and so that $\mathbf k\oplus \mathbf i$ could as well
have been written as $\mathbf k\ominus \mathbf i$.


Measuring noise floor in an impulse response Say I have an impulse response...

(an envelope of the actual IR obtained using the Hilbert transform)

... question: how do I calculate the surrounding noise level?
Visually this looks to be ~ -54 dB w.r.t. the peak.

But before you say something like "just average the surrounding parts" (ok maybe with better technical terminology), I would actually like to determine the noise level programmatically, in order to deal with wider peaks, like...

Just arbitrarily taking some noise samples could accidentally include part of the peak and obviously invalidate the result.

My hunch: set a threshold as the minimum number of samples above some level  X such that this level can be considered to be noisy. This could work well in the first case - as you increase X to above ~ -54 dB the number of samples above suddenly drops by a significant amount.

This approach may not work so well in the second case as the decay range consists of more samples (the increase may not be so significant). I could extend the sampling time... but I can only hold my breath for so long when recording.

I can suggest you take a look at the supporting information for Dirac software, where they describe the process of INR (IR to Noise Ratio) calculation.
The simplest approach is to estimate the noise level based on the beginning of the IR before the main peak. For that purpose, you can also do the tail of the Impulse response.

Yet another article supporting that estimation is by Lundeby et al.: "Uncertainties of Measurements in Room Acoustics". Their algorithm is an iterative procedure which tries to estimate the background noise level and the slope:

You could also try to verify that the noise floor doesn't have any distortions, by fitting the linear regression and making sure that the slope is as close to 0 as possible.


How can I detect when a signal deviates from a flatline? I apologize if this is not the correct forum, but I thought I'd ask.

I have an array of signal samples and I want to detect (retroactively, once I have measured the time/magnitude values) if there are any spikes / deviations during a sampled real signal
How would I do this intelligently?

I've been computing the mean of all my sets and the standard deviation and seeing if any of the individual samples greatly deviate from the mean, but I'm thinking there must be a more established method or approach

This type of problem is part of Detection Theory.  Do you have an estimate of the background noise levels?  How about the statistics of the noise?  If the noise is AWGN, things get a little easier.  You can set a detection threshold based off your noise statistics to achieve a desired probability of detection / probability of false alarm.

Be careful about using the mean of your data to set the threshold.  The "spikes" could bias your mean estimate.  If the noise is not stationary (i.e. it changes levels), using the mean could result in false detections (alarms).  One simple option is to compute a running average and create a threshold off this running average.  If the sample of interest is above the threshold, declare a detection and do not include the sample in the running average calculation.  There're a ton of variations to this simple concept depending on the situation.

If your signal is truly a flat line without noise, you could apply a highpass filter by subtract the previous sample from the current sample.  The result would be an array of zeros with non-zero sample values at the spikes.

For a more formal treatment of detection theory, I highly recommend Steven Kay's book "Fundamentals of Statistical Signal Processing, Volume II: Detection Theory".


Why do third-octave bands have a ratio of cube root of two? Without sounding tautological, I guess it boils down to the definition of the third octave, which is not: 2/3*Fc

but defined as the cube root of 2 ratio between frequencies.

Why is a third-octave not a "third of double the frequency" and instead "the cube root"? How does the square root play into this? Is there a log base 2 upon which this definition of a third octave is based?

To obtain the center frequency, it's the geometric mean of the upper and lower bounds, which gives you your addition sqrt 2, which gives you the upper bound Fc * 2^(1/6) = 2^(1/3) * sqrt(2) = 2^(1/3*1/2) and Fc / 2^(1/6) for the lower bound. Most tables start with 1kHz

A fixed interval corresponds to a fixed ratio of frequencies (at least in equal temperament, which is the most common tuning system in Western music). One octave is a factor of $2$. One semi-tone is a factor of $\sqrt[12]{2}$ because there are $12$ semi-tones in one octave. So one third of an octave (i.e., $4$ semi-tones) corresponds to a factor of $\sqrt[3]{2}$ because $\sqrt[3]{2}\cdot \sqrt[3]{2}\cdot \sqrt[3]{2}=2$ (i.e., one octave).


How can I learn image processing in a practical way? I am a beginner in image processing and I have to learn to program in MATLAB as soon as possible. I started with Rafael Gonzalez book, but understanding this theoretical knowledge and mathematical aspect from one side and bringing this knowledge into practice for doing my project from the other sides, it is a challenging task for me.I am working on 3D medical images. So I should learn how to do pattern recognition with these images.

Do you have any tips on how to translate my theoretical knowledge into practical knowledge?
Which programming language is suitable and convenient for learning image processing/computer vision and doing practical projects?

The book you mention in your question is a highly practical one. In fact since digital image processing (DIP) is by definition an applied field, any book on DIP must be and will be practical. They should be showing results of algorithms on images. But if you mean actual code, very little number of DIP books would give away anything other than pseudo or matlab / octave code.

In my opinion the preferred language for research in DIP is Matlab / Octave or Python, and the standard tool for production environment is C/C++ with MSVS (licence) or a number of excelent free alternatives.

Note that pattern recognition and medical image processing are two sub-fields of signal processing and image processing respectively. The former subject requires a lot of statistics and probability theory background, which is not included in a typical DIP book.


There is a mysterious factor when compute coeffiencts everyone
i m trying write higher order butterworth algorithm in c++.I can compute nominator and denominator second order. So i tell myself , i can cascade struct for higher order. But there is between factor two cascade and fourth order.

For example;

second order nominator a = [0,0976307608403291 0,195261521680658   0,0976307608403291] ,

denominator b = [1 -0,942806027702107   0,333329071063423]

when i do conv(a,a) in matlab result is

new_a  = [0,00953176546226155   0,0381270618490462  0,0571905927735693  0,0381270618490462  0,00953176546226155]

but i get this coeffs when i do butter func in matlab

butter_a = [0,0102094807912031    0,0408379231648125  0,0612568847472187  0,0408379231648125  0,0102094807912031]

maybe you see that between constant factor butter_a and new_a and it is 1.071 like that. I wonder what is that and how can i compute this variable for more than one cascading struct ?

sorry for my bad language

It sounds like you're comparing the coefficients for two different filters:

A fourth-order Butterworth filter designed for cutoff frequency $f_c$
A cascade of two second-order Butterworth filters, each designed for cutoff frequency $f_c$.

The two are not equivalent to one another, although in the passband and stopband you will see similar magnitude responses. Specifically, at the cutoff frequency $f_c$, you'll see 6 dB of attenuation from filter #1, while filter #2 would give you 12 dB. Therefore, it's not unexpected that you would get slightly different coefficients for the two implementations.


Memory of a filter used on real time digital signal processing My question is regarding the memory of a filter when we want to implement it on a real time dsp environment. We have certain memories like scratch memory , intermediate memory and internal memory. Does anyone knows how to we implement a filter in real time environment keeping in mind the memory considerations because memory is an important issue when it comes to real time signal processing. I want to implement a low pass filter using the above conditions.

If memory is your primary concern then an IIR filter is the way to go.  It uses very little memory and is pretty computationally efficient.

FIR filters' memory length is one less than the number of filter taps.  This usually isn't that much, but it can add up if there are a number of filters and/or they are long.


Different values when calculating signal's energy by hand vs MATLAB Given the following signal:
$$
x(n) =
\left\{
\begin{array}{ll}
      u(0.01n-0.025)2.3704e^{(-0.287682n)},   & n\ge 0 \\
      u(-(0.01n+0.025))2.3704e^{(0.287682n)}, & n< 0 \\
\end{array}
\right.
$$
I'm finding the energy using the summation formula:
$$E_{x(n)}=\sum_{n=-\infty}^\infty |x(n)|^2=\sum_{n=0}^\infty(2.3704e^{(-0.287682n)}u(0.01n-0.025))^2+\sum_{n=-\infty}^0 (2.3704e^{(0.287682n)}u(-(0.01n+0.025)))^2 =\sum_{n=0}^\infty(2.3704e^{(-0.287682n)}u(0.01n-0.025))^2+\sum_{n=0}^\infty(2.3704e^{(-0.287682n)}u(-(0.01n+0.025)))^2 = 2\times2.3704^2\times \sum_{n=0}^\infty (e^{(-0.287682\times 2)}u(0.01n-0.025))^n$$
So due to the sampling frequency and other information given about the problem, the sampling doesn't start until $n=3$ so I adjusted the limits of summation and removed the u(n) function:
$$=2\times2.3704^2\times \sum_{n=3}^\infty (e^{(-0.287682\times 2)})^n$$
And then adjusted the series to begin at $n=1$ so I could use the summation formula:
$$=2\times2.3704^2\times(\sum_{n=1}^\infty (e^{(-0.287682\times 2)})^n - (e^{(-0.287682\times 2)})^1 - (e^{(-0.287682\times 2)})^2)$$
So the series converges since the common ratio $<1$ and the sum of the series is:
$$S_{\infty}=\frac{a_1}{1-r}=\frac{(e^{(-0.287682\times 2)})^1}{1-e^{(-0.287682\times 2)}}=1.2857$$
Replacing the number above in the energy summation equation, I get:
$$E_{x(n)}=2\times 2.3704^2\times(1.2857-0.5625-0.3164)=4.5715$$
However, when I use MATLAB to find the energy of the signal using the following code:
E = 0
for i=1:fs+1
    E = E + abs(x(i))^2
end

where fs is the sampling frequency, the code yields the following result:
$$E=2.2858$$
which doesn't match my mathematical method. However, I noticed that the answer that MATLAB gives me is half the answer of what I got by hand. I also know with certainty that the MATLAB answer is correct.
So how do I go from $4.5715$ to $2.2858$? What's the step/reasoning that I'm missing?
Thank you!

Since I don’t know rest of the given code, I’ve decided to calculate the energy other way round.
The MATLAB code for calculating the given function’s energy is given below:
syms n; % Symbol assignment for summation

Energy_1 = symsum(power((2.3704 * exp(0.287682 * n)), 2), n, -Inf, -3); % Energy for the left-hand side of the sample domain

Energy_2 = symsum(power((2.3704 * exp(-0.287682 * n)), 2), n, 3, Inf ); % Energy for the right-hand side of the sample domain

Energy = Energy_1 + Energy_2 % Overall energy

According to this code above, the resultant overall energy becomes $E = 4.5715$.
As the number of samples increases, the result gets accurate. Also, the coefficient in the exponential’s power affects the accuracy of the result too, but this is not the main case.


Detecting changes between two similar images I have the following problem:

I need to detect changes between two more or less similar image.

In that case, I want to detect the missing parts.

Lets assume that the displacement is no more that 2px, and there is no significant change in light (but there is some) or camera pose.

What I want as an output, is a binary image about if there is a signifiacant change at that position.

My naive idea was, to apply blur (aroud the size of expected displacement), and compare pixelwisely, but in practice its not working really well. (mainly it has a huge amount of noise around edges)

Any better idea for this problem?

The first part of the problem that you're describing is called registration.
Assuming that the displacement is indeed small, and your problem has no time limits, you can use normalized cross-correlation on the images to find the displacement. It will also take care of amount of light normalization. (That means that you are assuming that the transformation between the two images is only displacement.)

After you found the displacement, you can normalize both of the images by their mean intensity, and do some simple pixel-by-pixel comparison.

Ideas for improvement and variations:

Use keypoints descriptors and matching to find transformation between two images, to save runtime.
Use different models of transformation for different scenarios - rotation for in-plane camera rotation that is orthogonal, scale for zooming scenarios, perspective for out-of-plane rotations.
Use smarter comparison between the two registered images, for example, only in the places where there are edges.
Create a noise model by checking your camera characteristics in order to know what is a significant change. For example, more than two standard deviations of noise.


Linearity, Causality and Stability of a System Consider a system:

$$
y[n] = y[n-1] + u[n],
$$

where $y[n]$ is the output and $u[n]$ is the unit step function.

Is this system causal, linear, time-invariant and stable ?

My attempt at the question is that for the systems without any input, we can never talk about all the properties above as they are related to input/output relationships of the system, whereas there are is no input in the system.

However, I'm inclined to think about the system as follows:

$$
y[n] = y[n-1] + \sum_{i=-\infty}^{\infty} 0 . x[i] + u[n].
$$

So trivially, the system is non-linear as it produces non-zero output for zero input, non-causal as it depends on future inputs, unstable as no matter the input, the output is unbounded and time-variant as if we delay the input, the output stays the same.

According to the interpratation (that belongs to you) that output $y[n]$ is given by the recursion relation $$y[n] = y[n-1] + u[n]$$ no matter what the input $x[n]$ is, then the following conclusions can be made about the linearity, time-invariance, stability and causality of the system.

First let's solve the output for $y[n]$

i.For $n<0$, since $u[n]=0$, we have $$y[n] = y[n-1]$$ which implies that $y[-1]=y[-2]=...y[-N]...= K$ for all $n < 0$

ii.For $n\geq0$, since $u[n]=1$, then we have
$$y[n] = y[n-1] + 1$$ which can be recursively solved to get

\begin{align}
y[0] &= y[-1] + 1 = K +1 \\
y[1] &= y[0] + 1 = (K+1) +1 = K+2 \\
y[2] &= y[1] + 1 = K +3 \\
&...\\
y[n] &= y[n-1] + 1 = K +(n+1) \\
\end{align}

Combining them yields: $$\boxed{y[n] = K + (n+1) u[n] }$$ for all $n$ as th fixed-unique output of the system for all inputs $x[n]$

Then we analyse the system for its properties:

The system is non-linear:
It's so evident since the system has a fixed, unique output for every input
it's a non-linear system as it cannot satisfy $$\mathcal{T}\{a x_1[n]+ b x_2[n] \} = a \mathcal{T}\{x_1[n] \} + b\mathcal{T}\{x_2[n]\} $$
The system is time-varying:
Again shifting the input have no effect on the output (which is fixed) so it cannot satisfy $$ y[n-d] = \mathcal{T}\{x[n-d] \}$$
NOTE: Does this make any sense? The system is already a fixed one. It has a fixed output for every input. Hence effectively nothing varies in its output. Nevertheless I called it as time-varying, because evidently it cannot saitsfy the time-invariance constraint.
The system is causal: it only depends on the past value of the output and it does not depend on any value of the input so it does not violate causality.
The system is unstable: as for a bounded input $x[n] = \delta[n]$ the output is $y[n] = K + (n+1) u[n]$ which grows unbounded as $n$ goes to infinity.
The system is not invertible.
The system can be implemented with and without any memory requirement.


How Does Color Conversion and Downsampling work When JPEG compressing a grayscale image? For color images I know you have to convert to YCbCr format, and then for downsampling is based on the fact that human beings see differences in Y more than differences in Cb or Cr.
How do I do these things but for grayscale (or monochrome) images?

Normally, monochrome images should be stored as images with a single channel. However, I have seen digital images with three channels containing the same values. This results in a displayed image looking grayscale.
Before compression multichannel data, it is common to decorrelate the channels. For RGB, one often concert them to luminance and chrominances. Several formulae can be seen, linear and nonlinear, to convert $R$, $G$ and $B$ to $L$ and the two $C_x$:
$$ L= \mathrm{c2l}(R,G,B)$$
and
$$ C_x= \mathrm{c2r}_x(R,G,B)$$
In most cases I know, for the same repeated channel (say $R$, but whatever),
$$ L= \mathrm{c2l}(R,R,R)\simeq R$$
and
$$ C_x= \mathrm{c2r}_x(R,R,R)\simeq 0$$
so the compression code will naturally "cancel" the unnecessary channel triplication for a grayscale image. As a side node, this can be the cause of errors in computing compression ratios,  as the original image is three times the size it should have been. This uses to happen  in  monaural (mono) audio stored on two duplicated stereo channels.
To be more concrete, with the classical linear YUV conversion:
Y =  0.299R + 0.587G + 0.114B
U = -0.147R - 0.289G + 0.436B
V =  0.615R - 0.515G - 0.100B

and since:
 0.299 + 0.587 + 0.114 = 1
-0.147 - 0.289 + 0.436 = 0
 0.615 - 0.515 - 0.100 = 0

one can see that the chrominance images are zero (or close to zero for nonlinear transforms), and remain very small in amplitude after DCT and subsampling, and will be unnoticeable in JPEG compression.


do i have an image processor to compare an image with a number of images and choose the identical image?? I require an image processor which compares an image with a set of images and finds the equivalent image out of them.

This is a very broad question and to be honest there are a multitude of ways to go about this. I highly recommend that you read more into imaging processing literature in general. From this you could perhaps suggest an approach that we could base our answers on, and help hone in on a good method for you.

That being said, a common approach to your problem is the SIFT algorithm. This algorithm will give back a set of features for a given image. From this you could construct a way to recognise images. However, this might not be the simplest solution if you are just a beginner at this game.

Best of luck!


Bilateral $\mathcal Z$-transform of exponential We all know that $a^nu(n)$ has unilateral $\mathcal Z$-transform. But what is the $\mathcal Z$-transform of $a^n$? (bilateral) When i tried to solve, i got answer as 'zero'.

But bilateral Laplace transform of $e^t$ doesn't exist. Both are exponentials  in discrete and continuous  domain respectively. Considering the similarity between Laplace and $\mathcal Z$-transform, how to explain the above problem?

Below, this is how I got 'zero'

$$a^n=a^nu(n) + a^nu(-n-1),$$

Now taking $\mathcal Z$-transform on both sides we get
$$\frac{z}{z-a}\quad \text{and}\quad \frac{-z}{z-a}$$ respectively which add to 'zero'

In complete analogy with the bilateral Laplace transform of $x(t)=e^{-at}$ (which doesn't exist), the bilateral $\mathcal{Z}$-transform of $a^n$ doesn't exist either. The series

$$\sum_{n=-\infty}^{\infty}a^nz^{-n}$$

converges nowhere, simply because $a^n$ grows without bounds for $n\rightarrow -\infty$ if $|a|<1$, or for $n\rightarrow\infty$ if $|a|>1$. Of course, for $|a|=1$ there series doesn't converge either.

EDIT:

As for your computation of the $\mathcal{Z}$-transform of $a^n$, the mistake lies in the fact that in addition to the algebraic expression of the transform you also need to consider the region of convergence. If you split $a^n$ (as you did) as

$$a^n=a^nu[n]+a^nu[-n-1]\tag{1}$$

you can compute the $\mathcal{Z}$-transform of both right-hand side expressions separately:

$$\mathcal{Z}\{a^nu[n]\}=\frac{z}{z-a},\quad |z|>|a|\\
\mathcal{Z}\{a^nu[-n-1]\}=-\frac{z}{z-a},\quad |z|<|a|\tag{2}$$

Note that the region of convergence (ROC) for the first part is outside the circle with radius $|a|$, whereas the ROC of the second part is inside the circle with radius $|a|$. The ROC of the total expression would be the overlap of the two ROCs, which is zero. Consequently, the sum doesn't converge anywhere and the $\mathcal{Z}$-transform of the total expression doesn't exist.


Autoregressive modeling (linear prediction) of electrical transmission lines? I've read that the reflection coefficients in speech processing (as computed by the Levinson-Durbin algorithm for solving the Yule-Walker equations) "represent the fraction of energy reflected back" at each tube junction,[1] assuming the speaker's vocal tract is modeled as a series of uniform lossless acoustic tubes (see Figure 1). Now, I've seen that statement from multiple PDFs online, but did not see it mentioned in my textbook (Rabiner, Schafer - Digital Processing of Speech Signals).

Does this mean we can do autoregressive modeling of electrical transmission lines as well? In other words, does a series of lossless transmission lines form an all-pole filter whose (magnitude) reflection coefficients (one for each junction) and linear-prediction coefficients can be estimated by inputting white noise to the system and then analyzing the output with the Levinson-Durbin algorithm? (Figure 2.)

Figure 1. Approximating a non-uniform acoustic tube using a series of short tubes of uniform cross-section. (Source: Rabiner, Schafer - Digital Processing of Speech Signals p.83.)

Figure 2. A series of lossless transmission lines. (Source: Pozar - Microwave Engineering 4th Ed. p.251.)

[1] http://research.cs.tamu.edu/prism/lectures/sp/l7.pdf

This is a neat idea and could probably be made to work - it seems like a variant on time-domain reflectometry, which exploits correlations of this kind https://en.wikipedia.org/wiki/Time-domain_reflectometry. However, there's an important difference between your two use cases: we can normally inject a test signal onto an electrical transmission line, while we normally can't do that in the vocal tract. Measuring the cross-correlation of the injected and reflected signals is generally going to be more reliable than measuring the cross-correlation of noise.


Fundamental misunderstanding of windowing I believe I have a fundamental misunderstanding of windowing. I thought that windowing is applied to a signal with pointwise multiplication, then some sort of filtering is then done to it (say using convolution in the time domain). However, looking back at my previous assignment, I seem to have applied the window to the filter itself(might be remembering incorrectly), which effectively means I convolved the signal with the window? Also, for images, it seems that I need to apply the window to a 2d kernel and not the image itself. This kernel then gets convolved with the image which means that I am once again convolving the image with the window, instead of pointwise multiplying. I understand that windowing is done to prevent spectral leakage/aliasing, but I dont quite get why it is being applied to the filter than the image itself.

If anyone could correct my misunderstanding I would be very grateful.

Your description of windowing a 2D image filter kernel, fits into a discrete-time FIR filter design method. Hence I will assume you are mixing windowing as a filter design method vs windowing applied to signals for Short-Time Fourier analysis or windowing applied for block-based processing of long signals. So please make it clear which action you are performing and the context of your intentions.

Then let me briefly describe windowing as a filter design tool.

An ideal discrete-time frequency selective LTI filter will have an infinite length and noncausal impulse response $h_i[n]$. So it shall be truncated and delayed for any practical FIR convolution applications.

Given an ideal filter impulse response $h_i[n]$, you obtain the coefficients of the practical impulse response $h[n]$ of length N samples via

$$h[n] = w[n] h_i[n]$$

Which is the point wise multiplication of the filter $h_i[n]$ with a window sequence $w[n]$ of length N. The main purposes of the window is to first truncate the infinite length ideal impulse response, and then to try to cure the side effects on the resulting filter spectrum caused by the crude truncation.

There are many different types of windows such as box-car,triangular, Hamming, Hanning, Blackman, Kaiser. The box-car is the crude rectangular window without any tapering. Other windows perform varying degrees of tapering at the ends of the impulse response so as to achieve a tradeoff between filter transition bandwith and passpand / stopband peak ripples.

Note that the windows's frequency response $W(e^{j\omega})$ is convolved with the ideal filter frequency response $H_i(e^{j\omega})$ which yields the frequency response of the practical filter;i.e,

$$H(e^{j\omega}) = \frac{1}{2\pi} H_i(e^{j\omega}) \star W(e^{j\omega}) $$.


First-Order Hold Filter Gain I want to know the amplitude of the first-order hold filter at the Nyquist frequency (the roll-off amplitude/gain).

I know that the Fourier transform of the reconstruction is given by:

$$\sum^\infty_{k=-\infty}F\left(\nu - \frac{k}{T}\right)\ \mathrm{sinc}^2(\nu T)$$

which corresponds to low-pass filtering of the spectrum in the frequency domain by a $\mathrm{sinc}^2$ function. This causes some aliasing due to the sinc's sidelobes. However at the Nyquist frequency $\nu = 1/2T$, approximately we have the amplitude:

$$\mathrm{sinc}^2\left(\frac{T}{2T}\right)\approx0.4053$$

Is this close to the correct value? Also, how should I display this gain as a dB reading? MATLAB's mag2db function gives -7.8, not sure if that's right.

That's correct, $0.4053$ is the approximate magnitude of the frequency response of linear interpolation at the Nyquist frequency. With the definition:

$$\operatorname{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$$

the exact value is:

$$\operatorname{sinc}^2\left(\frac{1}{2}\right)= \left(\frac{\sin\left(\frac{\pi}{2}\right)}{\frac{\pi}{2}}\right)^2 = \left(\frac{1}{\frac{\pi}{2}}\right)^2 = \left(\frac{2}{\pi}\right)^2 = \frac{4}{\pi^2}$$

or in dB:

$$20\log_{10}\left(\frac{4}{\pi^2}\right) = 20\frac{\log\left(\frac{4}{\pi^2}\right)}{\log(10)}\approx-7.844795081\text{ dB}$$


chirp signal - Frequency spectrum characteristics I am confused with my own intuition of the frequency response of a chirp signal.

In this case, the chirp signal is a frequency modulated signal with initial frequency starting at 100 Hz and ending at 300 Hz for 2 seconds.

Since the signal contains equal amounts of frequencies from 100Hz
to 300 Hz, the magnitudes of the frequencies should be equal? However, in the plot below, the frequencies at 100Hz and 300 Hz
have higher magnitude than the middle frequencies e.g. at 200Hz?
Should the plot be just a straight line?
Also, why does it appear that the middle frequencies are denser?
Is it correct to assume that since there are more frequencies in the
middle e.g. 100ish to 200ish than the end frequencies of 100Hz and
300Hz, then this should explain the increased density?

What is wrong with my assumptions?

Also do not neglect the fact that you are using a discretely sampled signal. You can create a 'cleaner' chirp by explicitly modifying the magnitude such that you have a perfectly clean magnitude and then convert it into the continuous signal using ifft. If you do fft of a signal generated with trig functions you will inevitably end up with artifacts due to the simple fact that your values are discretized and so may not fit exactly into the same fft bin when you transform them.


Image-processing - Dequantizer I begin to study how codec works, and I found a "summary-diagram" for DCT-codec (DCT stands for Discrete Cosine Transform). On this diagram, there are of course two steps: the encoding step and the decoding one. I try to make a reproduction of this diagram herebelow:

ENCODING STEP Image (8*8 blocks) -->| FDCT |-->| Quantizer |-->| Entropy Encoder |--> Compressed data

DECODING STEP Compressed data -->| Entropy Decoder |-->| Dequantizer |-->| IDCT |--> Image (8*8 blocks)

The diagram seems clear to me, except the "Dequantizer". I don't understand how the quantification step can be reversed, since it is this step that causes information loss. What does this "block" do? Can someone give me more precisions?

Thanks'.

Formally, the quantization is a mapping from a partitioned interval $[-A, A] = \cup_{k=1}^N I_k$ to the set of integers $\{1, \ldots, k\}$. The dequantization is the mapping from the quantization indices to a set of values $x_k$ chosen from each $I_k$.

When the quantizer used is uniform without any scaling (that is to say, when it is a simple rounding), the dequantizer is the identity. But when the quantizer is non-uniform, the integer passed into the entropy coder is not just a "rounded" value of the coefficient - an exponent, scaling factor, or any other monotonous transformation might have been applied to it. The dequantizer will have to reverse all the reversible steps in the quantification operation.

For example, in the AAC codec, the quantized value sent to the entropy coder is:

$$Q(x_i) = \text{sgn}(x) \left[\left(\frac{|x_i|}{s_i}\right)^\frac{3}{4} + 0.4054\right]$$

So the dequantizer will have to apply a $\frac{4}{3}$ exponent and multiply by the scale factor $s_i$ to recover the value.


LCCDE in simple words? What is LCCDE?I only know its abbreviation/full form :linear constant-coefficient difference equation

I know that in s domain we have differential equations and in z domain we have difference equation but what is LCCDE?What is its need?What is difference between a normal difference equation and a LCCDE?

An $N$th-order linear constant-coefficient difference equation (LCCDE) is of the form

$$y[n]=\sum_{k=0}^{M}b_kx[n-k]-\sum_{k=1}^{N}a_ky[n-k]\tag{1}$$

It is linear because the sequences $x[n]$ and $y[n]$ appear linearly in $(1)$, and it has constant coefficients because the coefficients $a_k$ and $b_k$ do not depend on the index $n$.

LCCDEs are important because they can be used to describe many practically useful discrete-time (or discrete-space) systems, such as linear time-invariant (LTI) filters (which are also called linear shift-invariant (LSI) filters, if $n$ does not represent a time index).

Given appropriate initial conditions, Eq. $(1)$ can be used to recursively compute the sequence $y[n]$ given its past values and the current and past values of the sequence $x[n]$. For causal LTI filters, $y[n]$ is interpreted as the filtered output, and $x[n]$ is the input. In that case, Eq. $(1)$ describes an infinite impulse response (IIR) filter.

If $a_k=0$, $k=1,\ldots,N$, Eq. $(1)$ reduces to

$$y[n]=\sum_{k=0}^{M}b_kx[n-k]\tag{2}$$

which is a convolution sum, and the coefficients $b_k$ represent the system's impulse response. Note that the impulse response has finite length $M+1$, i.e., the system described by $(2)$ is a finite impulse response (FIR) filter.


Modulation and filtering When I modulate a signal $x(t)$ with $\cos(2 \pi f t)$ and the modulated signal passes through a HPF, what output do i get in the frequency domain?

This is the core subject of linear modulation theory that's considered in analog communication systems. It's analysis is based on the modulation property of the Fourier trabsform.

Given a bandlimited baseband signal $x(t)$ such that its CTFT $X(\omega)$ is zero outside of the band  $-W < \omega < W$, its multiplication with the signal $\cos(\omega_0t)$ has the following spectrum:

$$ x(t) \cos(\omega_0 t) \longleftrightarrow \frac{1}{2\pi} X(\omega) \star \mathcal{F} \{ cos(\omega_0 t) \} $$ where $\star$ stands for convolution.

The Fourier transform of $\cos(\omega_0 t)$ is
$$  \cos(\omega_0 t) \longleftrightarrow \pi \delta(w-\omega_0) + \pi \delta(w+\omega_0)$$

Using the convolution property $ f(x) \star \delta(x-a) = f(x-a)$ i.e., signal shifts to the impulse position, we deduce the spectrum of the modulated signal as:

$$ x(t) \cos(\omega_0 t) \longleftrightarrow \frac{1}{2\pi}  \left( \pi X(\omega-\omega_0) + \pi X(\omega+\omega_0) \right) = 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0)$$

Now if you filter (convolve) this modulated signal with a (n ideal) BPF with a center frequency of $\omega_c = \omega_0$ and a bandwidth of at least $W_{BW} = 2 W$ then you would get the modulated signal spectrum alone. This can be shown by utilizing the filtering property of Fourier transform which says
$$ h(t) \star x(t) \longleftrightarrow H(\omega) X(\omega) $$

Now the defintion of the (ideal) BPF is such that its frequency response is:
$$
H(\omega) =
\begin{cases}
1 &, |\omega-\omega_0| < W \\
0 &, o.w. \\
\end{cases}
$$

Then the multiplication of $H(\omega)$ with $0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0)$ would yield:

$$ H(\omega) \left[ 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0) \right] = 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0) $$

Which means:
$$ \mathcal{F} \{ \text{BPF} \{ x(t) \cos(\omega_0 t) \} \}= 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0)$$

In this example the whole spectrum is the modulated spectrum so you don't actually filter out anything, but in other situations the BPF would indeed perform some useful operations.