Suggestions for creating 3D Bounding Boxes of Books on a Shelf from RGB-D Images I am new to image processing, so I beg your pardon if this is a trivial question. I am trying to use an RGB-D camera to detect 3D bounding boxes of books on a shelf.
Here are some examples of possible scenarios: books on a shelf 1 or books on a shelf 2. I would like to infer the dimensions (length, width and height) and the pose of each book and of the shelf. So approximate each of the element in the scene with a 3D bounding box.
I would like to get suggestions on which computer vision or image processing techniques are best suited to implement a solution to my problem. I am thinking of YOLO 5 or the following paper:

Deng, Z. and Jan Latecki, L., 2017. Amodal detection of 3d objects: Inferring 3d bounding boxes from 2d ones in rgb-depth images. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (pp. 5762-5770).

Are there any better or easier alternatives or already implemented packages?
Thank you all very much!

Exploiting the 3D geometry is important here. From the depth image, you could get the 3D point cloud by backprojection. Once the point cloud (set of 3D points) is given, you might benefit from VoteNet:

Deep Hough Voting for 3D Object Detection in Point Clouds, Charles R.
Qi, Or Litany, Kaiming He, Leonidas J. Guibas, ICCV 2019
https://arxiv.org/abs/1904.09664

If, in addition, the color (RGB) information is available, we could utilize within the same framework, this time maybe using ImVoteNet:

ImVoteNet: Boosting 3D Object Detection in Point Clouds with Image
Votes, Charles R. Qi, Xinlei Chen, Or Litany, Leonidas J.
Guibas, CVPR 2020
https://arxiv.org/abs/2001.10692

Both of the papers provide open source implementations. See here for VoteNet and this repository for ImVoteNet. Note that both of these are methods that require heavy supervision, i.e. large amounts of annotated training data. If such data is not available at your disposal, I would start looking into hand-crafted algorithms, e.g. descriptors or point pair feature detectors.
Hope this helps.


What are reasonable levels for the ERLE (Echo Return Loss Enhancement) in acoustic echo cancellation? I am working on acoustic echo cancellation which has a double talk detector. During the single talk, the ERLE (Echo Return Loss Enhancement) exceeds 100dB and during double talk around 20dB to 40 dB. I have read few papers and saw that for a good echo canceler the ERLE should be between 40-50dB. Am I going wrong anywhere? Is is ok to get such high ERLE?

It is annoying that I replied after 4 years, but this is helpful to others

Yes, ERLE has to be minimum of 40 dB for good echo cancellation. During the presence of only near end speech, you should not calculate ERLE as it will be near to 0dB, so ERLE has to be used only in presence of Far_end signal.

In your case I think you might be using low average filter length for calculation of ERLE. It mainly depends how you are calculating ERLE, ideally if you have ERLE 100 dB after convergence of your adaptive filter, then it is good enough.

Regards


Preamble length and channel estimation For OFDM we select the known gold sequence as the preamble or pilots. We can correctly estimate the channel(10 to 20taps).

Is there any relationship between the length of preamble and the number of taps of the channel?

For eg. Using a length of 8 pilot samples can we detect a channel of 20 taps?

Any insights?

For eg. Using a length of 8 pilot samples can we detect a channel of 20 taps?

If your OFDM system is sensibly designed, then, no, you need a pilot on every subcarrier. They don't necessarily have to be sent in the same OFDM symbol.

"Sensible" OFDM design means that you choose your subcarrier spacing no smaller than you have to. The reason that you're doing OFDM in the first place is always that you have a frequency-selective channel and would like to divide that big problem into many smaller problems by dividing that frequency-selective channel into ideally flat subcarriers. You'd call the bandwidth in which you can reliably consider the channel flat coherence bandwidth. If you're narrower than that, everything is fine, if your subcarrier is wider, you run into frequency selective fading trouble.

Now, you can't choose an arbitrarily large number of subcarriers – remember, that number is directly proportional to the temporal length of your OFDM symbol, and if you make that longer than it needs to be, you risk it becoming longer than the coherence time. That means your channel changes during your usage of it, and you run into temporal fading trouble.

So, these are competing constraints: you want as many subcarriers as possible for flatness reasons, but you want as little subcarrier as possible for constantness reasons.

Therefore, a well-designed OFDM system will have as little correlation between the channel coefficients of neighboring channels as possible, because if neighboring channels were reliably similar, then you'd be using more subcarriers than necessary.

Therefore, that system will simply have as many different subcarriers as it has subcarriers. You hence need a pilot on every single one of these, at least once in time.

Luckily, temporal changes tend to be relatively benign if your symbols aren't too long: phases tend to linearly shift, so extrapolation in time direction tends to work OK. Thus, modern OFDM systems tend to have "all pilot" symbols at the beginning of "superframes", and then pilots sprinkled across the time/frequency plane over the remaining OFDM symbols in that frame. That's sufficient to interpolate both in time as well as in frequency domain between these positions¹.

Look at LTE, for example (source of that image):

The red dots are the "sprinkled" pilots, and the green and yellow ones are "pilot-only"/synchronization OFDM symbols. The blue ones actually carry data.

¹ warning: he who interpolates OFDM symbols linear shall be disappointed by the performance compared to what he might have gotten using sinc interpolation.


Backprojection of 2D image on 3D model I am trying to implement the paper "Face Reconstruction in the Wild". In that I use a template 3D model and rotate it to align it with the orientation of the face in the 2D image. Now I want to project the 2D image onto the rotated 3D model. Could you please help me regarding how to do this. The paper gives no references for this so I think there might be an obvious solution I'm missing.

In section 2 of the paper it is mentioned that "Our image formation model assumes weak perspective
projection, relatively distant illumination, and Lambertian
reflectance" (emphasis is mine).

All three assumptions are making things much easier in various points of the algorithm further down the line but "weak perspective projection" makes the point of reprojecting the images on the model particularly easy.

A three dimensional projection is some $f: R^3 \rightarrow R^2$. There are essentially two types of it, orthographic and perspective. In an orthographic projection, points are projected vertically to the viewing plane. In a perspective transformation, there still is a viewing plane, but the actual rendered space shape is now the viewing frustum. Its shape (which maps to characteristics of a "camera") defines the final position of a 3D point on the 2D viewing plane.

A weak perspective transformation is essentially an orthographic but it involves a scaling too. What you scale with is the "average" distance of the object from the camera. So, given $x,y,z$ a point on some 3D mesh, its projection is roughly $P_x = \frac{x}{z}$. The bigger the $z$, the larger the scaling, the closer the point appears to the "ground line" (or horizon). When projecting "objects" (i.e collections of points linked by lines forming a mesh), you substitute the $z$ for the $z_{object}$ which is the average of all $z$ coordinates of all the points that make up the mesh. In this way, you still project vertically to the viewing plane but a far away object has a large $z_{object}$ which results in larger scaling of its $P_x, P_y$ coordinates which results in the object appearing smaller and closer to the horizon. Obviously, here, the object has to be somewhat coherent. A ball is ideal. A cube too. If you tried to render the star destroyer using this technique along other objects, then its perspective would look odd as its shape is such that it includes parts that extend very close to the camera but would have to be divided by the average $z_{object}$ of the whole mesh which is far, far away. Therefore, the "nose" of the destroyer would have to appear much bigger than what it will be rendered as.

This might seem like a digression but hopefully the spaceship example was illuminating on the characteristics of the weak perspective transformation because they are useful in pose matching.

Let's return to your case: The authors are assuming "weak perspective transformation". Essentially, they are saying that they assume that all portrait pictures have been acquired from far away enough to diminish the perspective effect of the camera lens because of the "relief" of the face (nose closer than eye socket). And that's great, because it means that you can use your fiducial points to pick up the orientation and scale of the face. You can then rotate (and scale) the 3D model and then project orthographically the image on to the 3D template. "Orthographically" as in, put the image on a plane, rotate it to the estimate pose around the 3D model and then extend lines vertically from the image / viewing plane towards the model. Once the line "touches" (intersects) the mesh, note down the $z$ position of that point (this now gives you the "depth" dimension you lost during taking the picture) and assign the image's pixel value to that point on the mesh (thus, projecting the image on the mesh).

Hope this helps.


Confused with filtering linear chirp signal using band-pass filter everyone. Recently, I am doing simulations about filtering a linearly chirp signal, and have been confused about the result. The signal was linearly chirped from 0 to 1000 Hz.
At first, the chirping rate is set to be 50 Hz/s. When the signal was filtered in time domain by a band-pass filter (Butterworth 500-600 Hz), the result was

--The first line of plot shows the signal

--the second line is the fft of the signal

--the third line is just the zoom in of the second line of plot.

We can see that the filtered signal corresponds to the pass band of the filter. 10 s corresponds to 500 Hz, and 12 corresponds to 600 Hz.

However, the the chirping rate is set to be 9850 Hz/s. the signal and the filtered signal was

since 500 ~ 600 Hz correspond to the time scale in 0.051 ~ 0.061 s (500/9850 ~ 600/9850). We can see the filter signal has been delayed and extended. Addition to that, we can see that the envelope of the filtered signal is oscillatory.

Is this normal or caused by improper parameter setup? Will I get the same result if the filter is an analog filter?

Any advice is greatly appreciated and thank you for your help.

You are seeing the difference between "wide band FM" and "narrow band FM" in your modulated waveform. Observe your modulation index for each case and then review the sideband levels versus mod index using Bessel functions. (The modulation index is the ratio of the frequency deviation over the frequency modulation). The first case is wide-band FM as your deviation (approximately +/-500 Hz) is significantly larger than your mod rate (50 Hz), so the mod index is >>1. The second case is narrow band FM as your deviation is significantly smaller than your mod rate (9850) and the mod index is <<1.  Narrow band FM will have only one pair of significant sidebands that will be at the mod rate offset from your carrier; while wideband FM will have a significant number of sidebands each spaced at your modulation rate and going out slightly beyond your deviation before becoming insignificant.

So when you are doing narrow band FM with one pair of sidebands, note that your signal will be dominantly your carrier with two sidebands that are each 9850 Hz away! You are seeing leakage products at your filter output; notice how small the amplitude is compared to the first case.  (Given your modulation is a ramp and not a sine wave, you will actually have multiple sidebands for each of the Fourier series components of your ramp signal, but the point is these will all be spaced by 9850 Hz from your carrier frequency which will also be dominant given the small moduation index).

While in your wideband case you will have spectral energy at spacing every 50 Hz with approximately 1KHz total bandwidth (review Carson's rule on this) and thus will provide tones that will pass through your filter.

Below is a plot of the Bessel function of the first kind from this site:
http://www.statsref.com/HTML/index.html?key_functions.html

This shows the relative level of the carrier (given by $J_0(x)$) and the first five sidebands $J_1(x), J_2(x),J_3(x),J_4(x), J_5(x) $ at a given modulation index x.


Can't understand multi intensity valued mask I am trying to do nuclei detection from a histopathological image. The original image and the mask (ground truth image or segmented image) are as follows:

The tissue image is (128,128,3) [RGB] and the mask is (128,128). Now, if I visualize with grayscale colormap:

My ultimate aim is to get bounding boxes for each nuclei in the mask.
What I have understood till this point: Do segmentation on tissue image, obtain the mask, find the bounding boxes on mask. Then I questioned why to get bounding boxes from mask, why not the original image? Then after some searching it led to conclusion that masks are more refined images only emphasizing required portions and setting remaining to 0.
Now, from my understanding a mask would be a binary image with certain pixels as 255 and rest as 0, so that when it is bitwise AND-ed with the original image we get only the required portions (wherever pixel values were 255).
So, in the given mask, how come there are different intensity values for nuclei? Is it due to the fact that some nuclei are strongly represented and some are weakly represented?

If so, then if I threshold the mask to 255 for every pixel value greater than 0, it should ideally become the mask which I described in my understanding, right?
If not so, then what is the meaning of this multi-intensity valued mask? It can't be applied like a binary mask for AND operations (since it has values like 38,60,110 apart from 255).
Even though it's a one channel image, the default colormap sets nuclei to different colors in the mask. Does it indicate that the mask is a m-ary mask with same color nuclei belonging to a certain class/label and if AND-ed with original image, then would yield only those nuclei of that class?

...masks are more refined images only emphasizing required portions and setting remaining to 0.

A mask is usually a binary bitmap that uses 1 at some location $(x,y)$ to indicate that it "selects" or "uses" or "includes" that pixel in the region of interest it defines.
The values that you see in your mask are not so much related to the mask but to its labeling.

...in the given mask, how come there are different intensity values for nuclei?

Each intensity value "tags" a different region. In other words, if your mask is in some $M$, then by evaluating $M == u$ where $u$ is a positive integer (say for example 5), then you would "select" all (connected) pixels that belong to the region of interest that has been labeled as "5".

Is it due to the fact that some nuclei are strongly represented and some are weakly represented?

No.
Also, "no" to #2 and "yes but not exactly" to #3: To get the behaviour described in #3 you would have to first "select" the mask (which casts it to a boolean) and then use it in logic expressions.
So, to get your bounding boxes all you have to do is iterate through the labels of the mask, extract it (with something like "where mask==label, set to 1 else 0" and find the min, max of the $x,y$ components of the pixel location.
Hope this helps.


What is the difference between an impulse response and a kernel? It seems to me they are two terms used to refer to the same concept, but, according to the Wikipedia page on Guassian filter

  a Gaussian filter is a filter whose impulse response is a Gaussian function

So, what is the difference between an impulse response and a kernel? I am particularly interested in knowing the difference in the context of computer vision and image processing.

The impulse response of an operation is the same as the kernel of that operation, but only if the operation can be represented by a kernel (the operation is linear and translation invariant).

For example, consider the following input (the impulse):

  [[0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 1., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.]]

If you give it as input to a laplace filter, you get:

  [[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]

This corresponds to the kernel of the laplace filter.

But if you give the same input to a 90 degrees rotation around the centre, you get:

  [[0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 1., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.]]

(the same as the input)

Obviously, this is not a 'rotation kernel', it is just the identity. This happens because there is no such thing as a rotation kernel, because rotation operations are not translation independent (besides the trivial case of rotating by 0 degrees).

If the operation was not linear, the kernel would not match the impulse response either (for example, for a thresholding operation).


What is the optimal Wiener-filter for this purpose? Suppose I want to measure something with a sensor which in response to $u(t)$ produces an output signal in form of $v_{out}(t) = (1-e^{-t/ \tau})u(t)$, where $u(t)$ is the Heaviside-function and $\tau=40ms$.

How to solve the deconvolution?

How do I determine the optimal Wiener-filter if I have an additional white noise $n(t)=n_0$ with $\omega_{max}=300s^{-1}$?

The noise which you defined will only shift your signal and only has zero frequency. i think you intended to define the noise's power. Also in addition to the noise power you need to know your signal's power spectrum to obtain Wiener filter.

Wiener filter will be (1/H)*(S/(S+N)). in this equation H is frequency response of your sensor, S is your signal's power spectrum and N is your noise power spectrum.

to obtain H you have to take derivative of your sensor unit step response to obtain unit impulse response, then take Fourier transform of unit impulse response.


On the meaning of s-plane and it's link to a transfer function Considering Fourier analysis and let's say I'm walking on the blue frequency axis in the below 3D plot from zero towards infinity:

So each time I encounter a non zero blue bar, I check the frequency at that point and measure the bar's height and say: Oh I found another sinusoidal component of the signal at this particular frequency and at this particular amplitude. This is how I grasp the idea of Fourier analysis.

But when I read about the s-plane I'm extremely confused to make it visual and link it to the actual signal as in the above fashion. I just cannot give a similar meaning what it is that really about.

Now imagine a system's output and input in time domain is given. Vout(t)/Vin(t) is transformed into s-domain and we call this ratio as H(s) which is the transfer function.

Here is a transfer function which is a surface on the s-plane:

Imagine I am walking on the s-plane. And imagine at a point on the s-plane I stop and the point s is a complex number and it is s = ωj + σ. Now if I look upwards towards the surface above me, obviously there will be a point on the surface. Lets call this point x. So what should it mean to me in terms of the signal or the transfer function?

Ok so at this particular point I found "what"? What is "x"? Is it a component of something? And my second question is if H(s) were a signal's transform not a transfer function what would x be in that case?

In the Laplace transform, being on the line $s = 0 \pm j\omega$ means you are looking at the pure sinusoids contained in the signal. This is exactly as you mentioned in your understanding of Fourier analysis. So this means that you are looking at the time-domain terms:

$$e^{j\omega{t}}$$

Which again represent the pure sinusoids composing your signal.

However once you are off this line and have some $s$ such that: $$s = \sigma \pm j\omega$$

Not only are you looking at the $e^{j\omega{t}}$ sinusoids terms, you are now looking at time domain terms:
$$e^{\sigma{t}}e^{j\omega{t}}$$

When talking about a stable system, poles being on the left-hand-plane of the $s$-domain means that $\sigma < 0$ so the sinusoids you are considering contained in the signal now have a damping envelope. If $\sigma$ was positive, the poles would be placed in the right-hand-plane and the signal would explode to infinity as time went on. This is why we want poles to stay out of the right-hand-plane if we want a stable system.

The takeaway is that in Fourier analysis you are considering the pure sinusoids composing your signal. With the Laplace transform you are considering sinusoids with that exponential term attached. In the case for a stable or oscillating system, you are considering the damped sinusoids ($s = -\sigma \pm j\omega$), decaying exponentials ($s = -\sigma$), pure sinusoids ($s = \pm j\omega$), and constant amplitude ($s = 0$) components composing your signal.

This is what makes Fourier analysis appropriate for identifying frequency content and the Laplace transform appropriate for analyzing stability and performance parameters.


Channel model: LOS Component and Rician K-factor I am asking about Rician K- Factor in this model:

I know that if $K=0$  then it is a pure Rayleigh fading, if $K= \infty$, then non fading Channel.

What is it, if $0<K<\infty$ ?

If i represet the NLOS component using the Kronecker model and take expected value from $\mathbf{H}$, then I get only loscomponent. Will be it right?

For $0 < K < \infty$, the channel is a combination of both a deterministic component (i.e., LOS) and a fading component.

As the $K$-factor is the ratio of the energy in the deterministic Line-of-Sight (LOS) component to the energy in the aggregation of the random scattered paths (i.e., the fading component), higher  means that the channel is more deterministic.


Derive Time-Domain Biquad Coefficients from Z-Domain I'm trying to derive the time-domain coefficients of a biquad filter based on its z-domain coefficients. The purpose is to import externally-generated EQ filters (in z-domain) into Max/MSP (in time-domain).
Could someone kindly explain how to do so?
Z-domain:
$$H(z) = \frac{b_0+b_1 z^{-1}+b_2 z^{-2}}{a_0+a_1 z^{-1}+a_2 z^{-2}}$$
Time-domain:
$$y[n] = a_0 x[n] + a_1 x[n-1] + a_2 x[n-2] - b_1 y[n-1] - b_2 y[n-2] $$

I'm not entirely sure what you are asking but if your Z-domain representation is
$$H(z) = \frac{b_0+b_1 z^{-1}+b_2 z^{-2}}{a_0+a_1 z^{-1}+a_2 z^{-2}}$$
Then your difference equation is
$$a_0 \cdot y[n] + a_1 \cdot y[n-1] + a_2 \cdot y[n-2] = b_0 \cdot x[n] + b_1 \cdot x[n-1] + b_2 \cdot x[n-2]$$
And solving for $y[n]$ we get
$$y[n] = \frac{1}{a_0}(b_0 \cdot x[n] + b_1 \cdot x[n-1] + b_2 \cdot x[n-2] - a_1 \cdot y[n-1] - a_2 \cdot y[n-2]) $$


How would I spectrally rotate/invert an audio signal (Matlab or Python)? So I have a 4.5 minute wav audio clip. I want to spectrally rotate the frequencies such that the spatio-temporal characteristics of natural speech are retained, while having the new audio be unintelligible. Any tips? I found https://ccrma.stanford.edu/~jos/sasp/Spectral_Rotation_Real_Signals.html, but I'm not getting anywhere with it.

$$ y[n] = (-1)^n \cdot x[n] $$
This will flip the spectrum in linear frequency. Easy to do, easy to recover.


demonstration using parseval I have to do a demonstration.
If we do the Parseval identity of the signals $x(t)$, $y(t)$ and $z(t)$ that go from $0$ to $T$ and that are real, we have:

$\int_{0}^{T} x(t)^2dt=\int_{-\infty}^{\infty}|X(f)|^2df$

$\int_{0}^{T} y(t)^2dt=\int_{-\infty}^{\infty}|Y(f)|^2df$

$\int_{0}^{T} z(t)^2dt=\int_{-\infty}^{\infty}|Z(f)|^2df$ (1)

If we sum up the these 3 equations, and we divide both sides by the interval $T$, we obtain:

$\frac{1}{T}\int_{0}^{T} (x(t)^2+y(t)^2+z(t)^2)dt=\frac{1}{T}\int_{-\infty}^{\infty}(|X(f)|^2+|Y(f)|^2+|Z(f)|^2)df$ (2)

Can you confirm me that everything is ok in this demonstration? I am trying to demonstrate this because my advisor says that since the PSD of the signals $x$, $y$, $z$, which are $S_{X}$, $S_{Y}$, $S_{Z}$, respect this relation:

$\int_{-\infty}^{\infty}(S_{X}(f)+S_{Y}(f)+S_{Z}(f))df=\frac{1}{T}\int_{0}^{T} (x(t)^2+y(t)^2+z(t)^2)dt$ (3)

then (and this is the point that is totally wrong for me) he says that $S_{X}+  S_{Y}+S_{Z}$ is the Fourier transform of  ( $x(t)^2 +y(t)^2+z(t)^2$):

$S_{X}+  S_{Y}+S_{Z} = \int_{-\infty}^{\infty}( x(t)^2 +y(t)^2+z(t)^2)\exp(-itf2\pi)dt$ (4)

this cannot be possible in my opinion. I already tried to show him the procedure that leads to the PSD of a signal etc. but he still believe (maybe because I am young, not expert and just arrived) that is true what he said.
So I am trying to do the above mentioned demonstration, saying that since $|X(f)|^2$ is different from the absolute value of the Fourier transform of $x(t)^2$, it cannot be contemporary true what it is stated in the eq. (4) and the application of the Parseval identity (2).

I beg you to tell me if my procedure is good. I want to convince him with something that is unassailable. (I am almost desperate, I have been trying to convince him that is wrong since 3 month ago)

Your first steps seem OK to me. Even if I'd put absolute values in the time domain as well (for generality), if the first three integrals exist, the fourth is valid. Integration is linear.

The weak point resides in the reverse inference. The equality of energy integrals does not mean that the integrands are related. For one, PDSs are real and positive functions. And there is NO reason why the Fourier transform of a sum of squared functions should be real and positive in general.

To convince your advisor, I'd suggest you to build your question on a clever counter-example.


Impulse Response to Frequency Response in Octave Having an impulse response of an audio system recorded as a wav file, how to calculate the frequency response of the system with octave?

The Frequency Response of a system is the Fourier Transform of the Impulse Response of the system.

Since we're in the real world and we have finite number of samples observed over finite time interval we use the Discrete Fourier Transform (DFT).
In Ocatve and MATLAB the DFT is implemented in efficient way using the Fast Fourier Transform (FFT).

In order to calculate it on MATLAB / Octave you can do the following:

% Assuming Impulse Response is given by 'vImpulseResponse'

vFreqResponse = fft(vImpulseResponse);

figure();
plot(abs(vFreqResponse));

Enjoy.


Difference between convolving before/after discretizing LTI systems Suppose I have transfer functions for two continuous causal linear-time invariant (LTI) systems:
$F_1(s)$ and $F_2(s)$.

Let $D\left\{\cdot\right\}$ denote the function that maps a transfer function from the continuous-time domain to the discrete-time domain (as another transfer function via the $z$-transform) using ZOH discretization.

What is the relationship between: $$D\left\{F_1(s) F_2(s)\right\} \quad \longleftrightarrow \quad D\left\{F_1(s)\right\}D\left\{F_2(s)\right\}$$

In other words, what is the difference between discretizing the convolution  of two transfer functions versus convolving their discretization?
Some brief experimentation has shown they are "almost" equal, but why?
Can we be precise about this difference?

As an aside,

Could someone point me to a reference on other properties of discretization, e.g., linearity $D\left\{a_1F_1(s) + a_2F_2(s)\right\} = a_1D\left\{F_1(s)\right\} + a_2D\left\{F_2(s)\right\}$?

Assuming you refer to the ZOH-discretization as shown in this figure from the mathworks site

the relationship between the continuous-time signals and the discrete-time signals can be derived as follows. The signal $u(t)$ is given by

$$u(t)=\sum_ku[k]g(t-kT)\tag{1}$$

where $g(t)$ is a rectangular impulse response (constant in the interval $t\in [0,T]$, zero anywhere else) representing the ZOH, and $T$ is the sampling period. Let $h(t)$ be the impulse response corresponding to the transfer function $H(s)$. Furthermore, let $f(t)$ be the convolution of $g(t)$ and $h(t)$:

$$f(t)=(g\star h)(t)\tag{2}$$

Then the continuous-time output signal $y(t)$ is given by

$$\begin{align}y(t)&=(u\star h)(t)\\
&=\int_{-\infty}^{\infty}u(\tau)h(t-\tau)d\tau\\&=\sum_ku[k]\int_{-\infty}^{\infty}g(\tau-kT)h(t-\tau)d\tau\\&=\sum_ku[k]f(t-kT)\tag{3}\end{align}$$

where I've used Equations $(1)$ and $(2)$. From $(3)$ the discrete-time output signal is easily obtained as

$$y[n]=y(nT)=\sum_ku[k]f((n-k)T)=\sum_ku[k]f[n-k]\tag{4}$$

with $f[n]=f(nT)$. So $y[n]$ is simply the discrete-time convolution of the input signal $u[n]$ and the total impulse response $f[n]$.

From this result it immediately follows that linearity must be satisfied, simply because convolution is a linear operation:

$$f(t)=\alpha_1(g\star h_1)(t)+\alpha_2(g\star h_2)(t)=(g\star (\alpha_1h_1+\alpha_2h_2))(t)\tag{5}$$

However, the (ZOH-)discretization of the concatenation of two transfer functions is not equivalent to the concatenation of the two discretized transfer functions. Referring to the figure above, the first case is equivalent to replacing $H(s)$ by the concatenation of $H_1(s)$ and $H_2(s)$. The other case involves concatenating two complete systems as shown in the figure, the first with $H(s)=H_1(s)$, the second with $H(s)=H_2(s)$. The difference between the two cases is that in the first case you only have one ZOH, whereas in the second one you get two ZOHs. The equivalent discrete-time impulse responses are

$$f[n]=(g\star h_1\star h_2)(nT)\tag{6}$$

and

$$f[n]=(g\star h_1)\star (g\star h_2)(nT)\tag{7}$$

which are generally not identical. The first case in Eq. $(6)$ is a ZOH discretization of the total transfer function $H(s)=H_1(s)H_2(s)$, whereas the second case in Eq. $(7)$ corresponds to a first-order hold (FOH) discretization of the total transfer function.


Basic Convolution Reverb Question I'm surveying for feasibility techniques for a potential future project that's rather out of my past areas of expertise, and I have a fairly basic question. Of course if it looks feasible I will have to learn a whole lot more on the topic.

What is the equation to compute the number of FFTs/IFFTs per sample required, given the number of samples in the FFT n, the number of samples in the FIR m, and the number of samples of audio data this is applied to at a time x?

Depends on your algorithm. Assuming overlap add/save and a filter length of N we would get:

Zero pad to 2*N and forward FFT: 4*N*log2(2*N) multiplies
Spectrum multiply: 8*N multiplies (4*N if input and filter are real)
Inverse FFT: same as 1
Overhead for padding, overlap handling, house keeping etc.

There are between 1 and 1.5 times the number of adds. There are also various ways to optimize this.

In general, convolution is not always the best way to create reverb. It is computationally expensive, it adds a lot of latency, and it is very difficult to update or modify in real time.


Sampling of a continuous function: Kronecker's or Dirac's delta? I've been reading some papers in signal proccesing and I'm very confused about the issue in the title of my question. Consider a continuous function of time $t$, $f(t)$, that I sample at uneven times $t_k$, where $k=1,2,...,N$. To me, it makes sense that the sampled function is:
$$f_s(t)=\sum_{k=1}^N\delta_{t,t_k}f(t),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)$$
where $\delta_{t,t_k}$ is Kronecker's delta (equals $1$ when $t=t_k$, zero elsewhere). However, in this paper, the author defines the sampled signal as:
$$f_s(t)=\frac{1}{N}\sum_{k=1}^{N}f(t)\delta(t-t_k),\ \ \ (2)$$
where $\delta(t-t_k)$ is Dirac's delta function and I really don't get why the $1/N$ appears here (the author claims that the sampling function is actually a weighted sum of delta functions
$$s(t)=C\frac{\sum_{k=1}^{N}w_k\delta(t-t_k)}{\sum_{k=1}^{N}w_k},$$
and here he choose $C=w_k=1$. I really didn't understand why). This last statement doesn't make much sense to me: the sampled signal would have infinite amplitude at $t=t_k$!

Despite of all this, it is much easier to define the Fourier Transform of $f_s(t)$ in the second case (equation $(2)$), because it is just the convolution of the window function (the FT of the Dirac comb) and the FT of the continuous signal $f(t)$, while on equation $(1)$ the FT is a bit more complicated because we have an integer function (Kronecker's delta) multiplied by a continuous function ($f(t)$). Any highlights on this?

Modeling of the sampling process via multiplication of a continuous-time signal by a train of Dirac impulses is the most common interpretation in my experience. If you dig into it deeply enough, you'll find some disagreement about the mathematical precision of this approach*, but I wouldn't worry about it; it's just a convenient model for the process. There are no impulse generators inside your cell phone's ADC generating periodic bolts of lightning that multiply their analog inputs.

As you noted, you can't calculate the continuous-time Fourier transform of the Kronecker delta function, as its domain isn't continuous (it is limited to the integers). The Dirac delta function, in contrast, has a simple Fourier transform, and the effect of multiplying a signal by a train of Dirac impulses is easy to show due to its sifting property.

*: As an example, if you're going to be mathematically precise, you would say that the Dirac delta isn't a function at all, but a distribution instead. But at an engineering level, these issues are really just semantics.

Edit: I'll address the comment below. You gave your mental model of the sampling process as:

$$
f_s(t) = \sum_{k=1}^N \int_{t_k-\epsilon_k}^{t_k+\epsilon_k} f(t) \delta(t-t_k) dt.
$$

The problem with this interpretation is that the typical ideal sampling model does not have that integration built in. Instead, it's a pure multiplication of the input signal by a Dirac impulse train. If you look more closely at the equation you showed for $f_s(t)$, you'll see that the right hand side actually has no independent variable; $t$ is the dummy variable of integration. For any $\epsilon_k > 0$ in the above, according to the Dirac impulse's sifting property, you would get:

$$
f_s(t)=\sum_{k=1}^N f(t_k),
$$

which isn't correct. Instead, the model for the sampled signal is:

$$
f_s(t)=\sum_{k=-\infty}^{\infty} f(t) \delta(t-kT)
$$

Which is very similar to the above, except generalizing for an infinitely long impulse train along the time axis and assuming that the data is uniformly sampled at time instants $t_k = kT$. The Fourier transform of the resulting signal is:

$$
\begin{align}
F_s(\omega) &= \int_{-\infty}^{\infty}f_s(t) e^{-j \omega t} dt \\
&= \int_{-\infty}^{\infty} \sum_{k=-\infty}^{\infty} f(t) \delta(t-kT) e^{-j \omega t} dt \\
&= \sum_{k=-\infty}^{\infty} \int_{-\infty}^{\infty} f(t) \delta(t-kT) e^{-j \omega t} dt \\
&= \sum_{k=-\infty}^{\infty} f(kT) e^{-j \omega kT}
\end{align}
$$

If we define the discrete, sampled version of the signal $f(t)$ to be $x[n] = f(nT)$, then you are left with:

$$
F_s(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n}
$$

which is exactly the definition of the discrete-time Fourier transform.


Information Loss in a Multirate System Assume $M=3$ and $L=2$ for the multirate system shown below. Furthermore, assume an arbitrary fourier transform $X(e^{j\Omega})$. This one lies between $-\pi$ and $\pi$. The fourier transform after the filter $H_{LP}(z)$, which is an ideal low-pass filter with a cutoff-frequency of $\Omega=\frac{\pi}{3}$, lies between $-\frac{\pi}{3}$ and $\frac{\pi}{3}$, since we sample up by a factor of three. If I then sample down by two, the resulting fourier transform lies between $-\frac{2\pi}{3}$ and $\frac{2\pi}{3}$. Now the spectra of the input and output will clearly look different, therefore also the signal in time. Is it right to say that with this specific multirate system there will be information loss, although no aliasing effects occur?

Note: In case my problem description is confusing, I can provide an example.

Since in your system there is no aliasing and you do not filter away important parts of the signal (because LP is having the correct cutoff-frequency), there is no information loss.

When a system changes the signal shape, it does not necessarily mean that information is lost. As long, as the operation is reversible, the same amount of information is contained in $y[l]$ and $x[n]$. However, even though the same information is contained, it can be differently difficult, how to extract the useful information from a given signal, depending on its shape.

However, in reality this ideal low-pass filter does not exist, and it might lose some information.


White Road Sign Detection Possible Duplicate:
  How to extract traffic signs from a photograph?

I am working on a project which involves detecting north american speed signs (black/white) using a 720p HD camera. My problem is that I cannot seem to figure out how to classify the rectangular profile of the speed signs in order to pass it on to the text detection stage. My current segmentation approach detects differing levels of white in an image and thresholds the image based on those values (determined from color of road signs). The problem is that I need to segment it further as threshold-ed shape doesn't contain the complete road sign, but bits and pieces of it. I wanted to know what would be the best way to classify the road sign?

These are some pre-segmented images and the corresponding segmented images:

As you can see the road sign has holes in it but its distinct rectangular profile is in-tact in most cases. I am currently looking into using kurtosis shape profile to determine if it is rectangle however the approach seems complicated and computationally expensive, I wanted to know if anyone knows a faster simpler method to determine the rectangular profile of the road sign from the segmented image?

I did this in Matlab using the image processing toolbox.

EDIT: Use Local OTSU thresholding it will allow for segmentation and further morph operations allow to speed sign detection.

Whoever marked this similar to the link below as a duplicate, its not... there is fundamental difference between chromatic and achromatic colors (This question deals with achromatic) whereas the link suggested below deals with chromatic.

Your question is broadly answered here: How to extract traffic signs from a photograph?

Your specific case is a little more complicated because the signs are not particularly distinguishing, however the same suggested techniques of cross correlating the images with clean images of original signs (you will probably want a library of images containing all signs for best results) are still applicable.

One suggestion above was to use an infrared camera which will be very effective if the road signs are micro prismatic or glass bead designs.


Quantization step-size for a given signal based on its PDF I am familiar with the principles of midtread and midrise quantizer. However, I have difficulties determining the step size where it hasn't been explicitly given. For example, the following probability distribution of the input signal s(t) is given. This signal shall be quantized with 2 bit and we want to use a uniform, midrise quantizer. How can I calculate the step size and thus reconstruction values/decision thresholds for this quantizer?

You have a uniform input PDF, and the optimal quantizer for a uniform inpu will be a uniform quantizer. Note that the your uniform PDF has two pieces and the interval from $-4$ to $4$ is completely untreated.

Then, assuming a 2-bit (4 level) uniform midrise quantizer with the given PDF your decision intervals $I_k$ and reconstruction levels $y_k$ will be as follows:

$$ x \in I_1 = [-6,-5] \implies y_1 = -5.5 $$
$$ x \in I_2 = [-5,-4] \implies y_2 = -4.5 $$
$$ x \in I_3 = [4,5] \implies y_3 = 4.5 $$
$$ x \in I_4 = [5,6] \implies y_4 = 5.5 $$


Create a mask like DAVIS dataset I am testing the dataset DAVIS from here with one of the papers it has. The codes works fine with DAVIS dataset but it does not work for my data. In particular, the code does not accept the mask I create for the images. Does anyone know how to create a mask like in the DAVIS dataset. Is there a recommended way for that? Can someone please help?

Here are a few tools you might find helpful for annotating and generating Ground Truth segmentation masks:

Xyonix Vannot - Video Annotation Tool for Object Segmentation
Best image annotation platforms for computer vision (+ an honest review of each)
Playment.io's Video Annotation tool


2-D sampling and reciprocal lattice I got this problem wrong, but I'm not sure why. I know that the FT of a signal is fully determined within 1 unit cell of a reciprocal lattice and outside this range are replicas. Here's the problem:

Edit: I didn't add it originally because I thought it wouldn't display well, but here's a screenshot of the exact solution I was given:

I dont know, if my attempt is correct, but here could be one argumentation (which actually conflicts with both solutions):

Intuitively, I would have agreed to your solution, because the spectrum is equally wide in height and width. But, then what about this argumentation:

The circles have radius $a/2$, and their center is at $\pm (\frac{a}{2},\frac{a}{2})$.
If I rotate the spectrum by 45 degrees clockwise, I get the center of the circles at $(\pm\frac{\sqrt{2}}{2}a,0)$. Then, the spectrum lies between $\pm\frac{1+\sqrt{2}}{2}a$ in x-dimension and $\pm \frac{a}{2}$ in y-direction.
In this case, I can use two different sampling frequencies in the new x- and y-direction, namely $a$ in the y-direction and $a(1+\sqrt{2})$ in the x-direction.
This idea could explain the non-symmetric sampling that your prof suggests. However, I do not get the same values as he/she does.


why multiply grayscale images by 256 I am a newbie in image processing. in an open source code for convolution neural networks, the author is multiplying a gray scale image from MNIST dataset by 256. Please  can any one tell me why?

ADDITION: the greyscale images have a max value of 0.8681 and a min of -0.1319

Many thanks.

Grayscale images are typically stored as 8 bits per pixel in files. Perhaps the network was trained to work with such data. $2^8 = 256,$ so 8 bits can represent any number from 0 to 255. It may well be that the author was converting from range $[0, 1]$ to $[0, 255]$. Note that $0.8681 - (-0.1319) = 1.$ Multiplying by 255 to stay within the range would have been more appropriate, but multiplying by 256 is appealing because it does not introduce any rounding error in finite precision floating point arithmetic.

Subtracting the mean does not change the story as it can be equivalently done before or after the multiplication.


How do bode plots work with unstable systems work? If I had a system with right-half s-plane poles, how would a frequency response work? Since a purely imaginary value for s, would cause the Laplace transform to diverge for such a system, what meaning would a bode plot have in such a situation?

The frequency response can exist, even if there are poles in the right half-plane, as long as there are no poles on the $j\omega$-axis. You cannot use the unilateral Laplace transform, instead you must use the bilateral Laplace transform (with the lower integration limit at $-\infty$), or, simply, the Fourier transform.

The corresponding time domain function is left-sided (or two-sided, if there are also poles in the left half-plane), i.e., non-causal. Note that poles in the right half-plane only indicate instability if you restrict yourself to causal systems. Non-causal systems can be stable if there are poles in the right half-plane.

Take as an example two simple first order systems, one with a pole in the left half-plane, and one with a pole in the right half-plane, such that it is a mirror image of the left half-plane pole:

$$H_1(s)=\frac{1}{s+1}\tag{1}$$
$$H_2(s)=\frac{1}{s-1}\tag{2}$$

If $H_1(s)$ is the frequency response of a causal system (note that we can decide that in this case), then it is stable, and its frequency response is simply

$$H_1(j\omega)=\frac{1}{j\omega+1}\tag{3}$$

In the same way, if we decide that $H_2(s)$ describes a non-causal but stable system, then its frequency response exists and is given by

$$H_2(j\omega)=\frac{1}{j\omega-1}\tag{4}$$

Note that their Bode magnitude plots are identical, because

$$\big|H_1(j\omega)\big|=\big|H_2(j\omega)\big|=\frac{1}{\sqrt{1+\omega^2}}\tag{5}$$

Their phase plots are of course different.


What is the math behind median filter's noise reduction property? I am interested in understanding the mathematical reason for why does applying a median filter on an image (or signal) result in reduction of noise.

Intuition: The intuition is this: Your noise is some event or events that are rare, and that when compared to other events, look like outliers that shouldn't really be there.

For example, if you are measuring the speeds of every car on the highway as they pass by you and plot them, you will see that they are usually in the range of say, $50$ mph to $70$ mph. However as you are inspecting your data for your boss, you see that you recorded a speed of $1,000,000$ mph. Not only does this value not make physical sense for the speed of an actual car on a highway, but it also sticks out wildly from the rest of your data. Chalking this event up to some strange measurement error, you remove it, and give the rest of your data to your boss.

However as you continue your measurements day in and day out, you notice that every now and then, you get those wild measurements of speed. For example, over the span of 1 hour, you measure 1000 cars, and their speeds are nicely between $50$ and $70$ mph, however 3 of those have speeds of $23,424$ mph, $12,000,121$ mph, and $192,212,121,329,982,321,912$ mph, breaking not only local state laws, but also those of theoretical physics.

You get tired of continuously having to go in, and remove those errant data points caused by your cheapo-radar by hand. Afterall, your boss is really only interested in the statistics of the speeds, not so much every actual value. He likes to make nice histograms for his bosses.

Those errant and large numbers are a kind of 'noise' you reckon - 'noise' caused by your cheapo-radar that you bought from a shady pawn shop. Is the noise additive white gaussian noise? (AWGN). Yes and no - It's spectrum is wideband and white, but it is temporally rare, sparse, and very localized. It is better referred to as 'salt and pepper' noise, (especially in the image processing domain).

Thus, what you can do, is run your data through a median filter. Your median filter will take a block of say, $5$ speed points, (points 1 to 5), find the median, and spit that value as the 'average' speed. Then it will take the next 5 points, (points 2 to 6), take that median, and spit out this as the average, etc etc.

What happens when you come across one of your faster-than-light speeds?. Let us say that your 5 speeds were [45, 65, 50, 999999, 75]. If you took the normal average, your 'average' speed here will be something quite large. However if you take the median, your 'average' will be 65. Which best approximates the average that you are really trying to measure? The median metric.

Thus, if you filter your data with a median filter, you will be sure to remove those outliers - and you have thus faithfully 'de-noised' your signal. In contrast, if you tried to remove your noise via traditional filtering, (nothing but a moving weighted sum), you will instead 'smear' the error across your data, and not get rid of it.

Math: The math is this: The median measurement is what is referred to as an order statistic. That is, it returns the value of your data, along some point, after it has been ordered. The max and min are also both order statistics - they return the extreme points of your data after it has been ordered. Taking the median also returns the value of your ordered data, but right from the middle.

But why are they different from mean filters? Well, mean-filters compute an average using all the data. If you notice from max, min, and median, you are getting an answer without using all the data. In fact, all the median does is order your data, and pick the value in the middle. It never 'touches' the outliers, like those large speeds that you measured.

This is why median - an order statistic - is able to 'remove' outlier noise for you. Outlier noise segregates itself in front of the median, and the median never comes near it or considers it, while still giving you a nice estimate of central tendency.


Why LTI system cannot generate new frequencies? Why $Y (\omega) = X(\omega)H(\omega)$ implies that an LTI system cannot
generate any new frequencies?
Why if a system generates new frequencies, then it is not LTI?

One of the definitive features of LTI systems is that they cannot generate any new frequencies which are not already present in their inputs.
One way to see why this is so, comes by observing the output's Fourier transform $Y(\omega) = H(\omega)X(\omega)$
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau \longleftrightarrow Y(\omega)=X(\omega)H(\omega),$$ holds only when the impulse response $h(t)$  exists and it will exist only when the system is LTI.
From a little thought, guided by a simple graphical plot, and using the convolution-multiplication property above, one can see that the frequency region of support $R_y$ (set of frequencies for which $Y(\omega)$ is non-zero), of the output $Y(\omega)$ is given by the intersection of the regions of support $R_x$ and $R_h$ of the inputs $X(\omega)$ and frequency response $H(\omega)$ of the LTI system:
$$R_y = R_x \cap R_h$$
And from set algebra we know that if $A = B \cap C$ then $A \subset B$ and $A \subset C$ . That is, an intersection is always less or equivalent to what are being intersected. Therefore, the region of support for $Y(\omega)$ will be less than or at most equal to the support of $X(\omega)$. Hence no new frequencies will be observed at the output.
Since this property is a necessary condition for being an LTI system, any system that fails to posses it, therefore, cannot be LTI.


obtaning discrete-time system impulse response $h[n]$ from step response $s[n]$ The discrete-time system of figure below a step response: $$s[n] = \left[ 5 - 4 \left( \frac{4}{5} \right)^n  \right] u[n]$$

Find the impulse response of the system and verify it from the system
difference equation.
sr i want ans of this question.with explanation

Considering an LTI system, its impulse response $h[n]$ and the step response $s[n]$ are related by $$h[n] = s[n]-s[n-1]$$ and $$s[n] = u[n] \star h[n] = \sum_{k=-\infty}^{n} h[k] = \sum_{k=0}^{\infty} h[n-k] $$

Given your definition of $s[n]$ therefore its $h[n]$ is
$$h[n] = \left[5 - 4 \left( \frac{4}{5} \right)^n \right] u[n] -\left[5 - 4 \left( \frac{4}{5} \right)^{n-1} \right] u[n-1]$$

After simplifications it would become:
$$h[n] = 5\delta[n] + 5 \left( \frac{4}{5} \right)^{n} u[n-1] - 4 \left( \frac{4}{5} \right)^{n} u[n]$$ further simplifies to

$$ \boxed{ h[n] = 5 \left( \frac{4}{5} \right)^{n} u[n] - 4 \left( \frac{4}{5} \right)^{n} u[n] = \left( \frac{4}{5} \right)^{n} u[n] }$$

Looking at the signal flow graph and using Z-transform variables you may achieve the same result by considering:

$$Y(z) = \frac{4}{5} Y(z) z^{-1} + X(z)$$
$$Y(z)( 1 - \frac{4}{5} z^{-1}) = X(z)$$
$$ H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1 - \frac{4}{5} z^{-1}} $$

From which the inverse Z-transform (with the causal ROC) yields:
$$ \boxed{ h[n] = \left( \frac{4}{5} \right)^n u[n] }$$ which is the same result obtained from the time-domain approach.


OFDM Transmitter knowledge I have some questions regarding some concepts used in designing the OFDM transmitter which i got from different literature reviews but i got mixed up with these concepts.

•   The N-Ifft points and number of subcarriers are not same all the time. Why we introduce unmodulated or zero padded subcarrier between the modulated or useful subcarrier.

•   Should the training sequence be added before ifft operation or after ifft+cp inclusion operation? If both can be possible then what is the best approach between these two.

•   For the DAC conversion, raised cosine pulse shaping filters are considered to be the best among other shaping filters. What should be its roll off factors, No. of symbols and No of samples in each symbol (for OFDM). How much should be an oversampling factor for upsampling process.

•   Is Low Pass Filter should be used at the transmitter side too what is the purpose of using it at the transmitter side ?

I hope to get a prompt reply.

Thanks

• The N-Ifft points and number of subcarriers are not same all the time. Why we introduce unmodulated or zero padded subcarrier between the modulated or useful subcarrier.

For different reasons. Mainly:

Most OFDM systems are based on direct-conversion receivers. In those receivers, you mix a signal down to complex baseband by multiplying it with a complex tone of the signal's center frequency – but that leads to the signal at exactly 0Hz in the baseband to be unrecoverable, technologically. For that reasons, the DC carrier (or carriers around DC) (DC=direct current, a very Electrical Engineering way of saying "0Hz") is not used
OFDM systems often have to deal with problems of strong Out-Of-Band radiation. By leaving the "edge" carriers unused, you can mitigate that.

  • Should the training sequence be added before ifft operation or after ifft+cp inclusion operation? If both can be possible then what is the best approach between these two.

I assume you want to train your equalizer? There's a lot of approaches to channel equalization with OFDM, and every one has its advantages and disadvantages.

The same applies to adding a timing and frequency synchronization sequence: there's a lot of approaches to doing synchronization for OFDM, and hence, as many trade-offs for advantages and disadvantages to having something in time- and frequency domain. Note that fine frequency estimation is often done with known time-domain signal techniques (because we can't do an FFT overly well if we don't even know where), and coarse frequency correction (ie. correcting an offset of integer multiples of the subcarrier spacing) can then be done afterwards (in time- or frequency domain, whichever is more appropriate for the system.

I personally would recommend Schmidl&Cox' classic paper Robust frequency and timing synchronization for OFDM – it's a great read, and I personally think it's a good paper, with a short, concise intro, a well-explained method, well-condensed mathematical proofs, and useful explanation and derivation of performance parameters. Also, if you're doing OFDM in your thesis, and want to actually have read at least one of the papers you cite, this should probably be the one.

  • For the DAC conversion, raised cosine pulse shaping filters are considered to be the best among other shaping filters. What should be its roll off factors, No. of symbols and No of samples in each symbol (for OFDM). How much should be an oversampling factor for upsampling process.

Generally, as you've learned, Roll-off etc are parameters of your design that you choose to make your design work according to your requirements. So, no matter for which method you'd be asking, "what are the values of the Roll off factor for ?", the answer will always be, "what you design them to be!".

However: OFDM was partially "invented" to eliminate the need for a set of pulse shaping filters. So your question regarding a raised cosine doesn't make much sense. If you pulse shape the time domain OFDM signal, you destroy the "O" in OFDM! You mustn't do that.

  • Is Low Pass Filter should be used at the transmitter side too what is the purpose of using it at the transmitter side ?

I think you're sitting in front of some specific example, and think it is "all possible OFDM transceivers". I can't tell you what your example is – I don't know it. Maybe that low-pass filter is just part of the direct conversion mixer that converts the OFDM baseband signal to a passband signal?!


Impulse response of ideal filters I am aware that an ideal low-pass filter in both continuous time and discrete time has a $\mathrm{sinc}$ impulse response. What would the impulse response of an ideal high-pass or band-pass filter look like (assuming that it is real)? My intuition says that it would be some sort of superposition of shifted $\mathrm{sinc}$ functions but I am not quite sure.
Thanks in advance.

EDIT:
Here is my attempt at deriving the impulse response $h(t)$ for a continuous-time ideal high pass filter with a zero phase characteristic.

The frequency response in this case can be written as $H(\omega) = u(\omega - \omega_c) + u(-\omega -\omega_c)$ where $u$ is the unit step function and $\omega_c$ is the filter cut-off frequency.

Using the Fourier transform pair:

$$u(t) \leftrightarrow \frac{1}{j\omega} + \pi \delta(\omega)$$

we can exploit the duality property to obtain:

$$\frac{1}{jt} + \pi \delta(t) \leftrightarrow 2\pi u(-\omega)$$

which implies (from linearity):

$$\frac{1}{j2\pi t} + \frac{1}{2} \delta(t) \leftrightarrow u(-\omega)$$

Additionally we obtain from the time reversal property and the scaling property of the delta function:

$$-\frac{1}{j2\pi t} + \frac{1}{2} \delta(t) \leftrightarrow u(\omega)$$

Applying the inverse Fourier transform to:

$$H(\omega) = u(\omega - \omega_c) + u(-\omega -\omega_c) = u(\omega - \omega_c) + u(-(\omega +\omega_c))$$

and using the linearity property, the frequency shift property, and the two Fourier transform pairs I wrote above we obtain:

$$h(t) = e^{j\omega_c t} \{ -\frac{1}{j2\pi t} + \frac{1}{2} \delta(t) \} +  e^{-j\omega_c t} \{ \frac{1}{j2\pi t} + \frac{1}{2} \delta(t) \} $$

and I don't know where to go from there (specifically, I don't know how to handle the delta functions that popped up).

Given that you are ok with the impulse response of the ideal zero phase lowpass real filter, then you can easily derive the impulse response of the ideal high pass zero phase real filter. In continuous-time let $$h_{lpf}(t) = \frac{\sin(\omega_0 t) }{\pi t} $$ be the impulse reponse of an LTI ideal zero phase lowpass filter whose frequency response is by definition:
$$ H_{lpf}(\omega) = \begin{cases}  1 ,& \text{ for } |\omega| < \omega_0 \\ 0 ,& \text{ otherwise } \\ \end{cases} $$

Then you can define the frequency response of the ideal zero phase highpass filter as:
$$ H_{hpf}(\omega) = \begin{cases}  1 ,& \text{ for } |\omega| > \omega_0 \\ 0 ,& \text{ otherwise } \\ \end{cases} = 1 - H_{lpf}(\omega)$$

By taking the inverse Fourier transform of $H_{hpf}(\omega)$ you can conclude that
$$ \boxed{ h_{hpf}(t) = \mathcal{F}^{-1}\{ 1 - H_{lpf}(\omega) \} = \delta(t) - h_{lpf}(t) }$$

Similarly in discrete-time you can conclude that:
$$ \boxed{ h_{hpf}[n] = \mathcal{F}^{-1}\{ 1 - H_{lpf}(e^{j\omega}) \} = \delta[n] - h_{lpf}[n] }$$


What is the difference between Kalman filter algorithm and stationary Kalman filter algorithm? I want to compute the stationary Kalman filter algorithm but I haven't found any information about that algorithm ( not even the pseudo code ) so, I wonder what is the difference between the Kalman filter algorithm and the stationary Kalman filter algorithm and how I can derive it from the Kalman filter algorithm. Could someone point me out resources I can read about it?

I just found this post:

How to derive the stationary Kalman filter predictor?

and this document:

http://www.uh.edu/~bsorense/kalman.pdf on page 5.

Basically what I understand is that you calculate beforehand ( before the recursive process ) the Kalman gain and the state covariance matrix when the Kalman gain goes to infinite using these two equations:

$$\bar{P} = A\bar{P}A^T - A\bar{P}C^T ( C\bar{P}C^T + R )^{-1}C\bar{P}A^T + Q$$
$$\bar{K} = \bar{P}C^T(C\bar{P}C^T+R)^{-1}$$

and in the recursive function, this equation:

$$\hat{x}(t+1|t) = (A-A\bar{K}C)\hat{x}(t|t-1) + A\bar{K}y(t)$$

Is that correct?

Notes:

Yes, it is a question from a DSP book ( https://users.aalto.fi/~ssarkka/pub/cup_book_online_20131111.pdf excercise 4.6 on page 83 ) but there is no information in that book about the stationary Kalman filter algorithm.

I’m not sure what you mean by the stationary Kalman filter, but it seems to be what I would call the steady-state Kalman filter.

If that is the same thing, then you just solve for the Kalman gain at $t=\infty$ and apply the normal Kalman filter equations.

I’m not sure you mean “the Kalman gain goes to infinity” in a well-posed problem.


Whats is the difference between FIR/IIR filters and Chebyshev/Butterworth filters I am new to Signal Processing. From my understanding -- FIR/IIR just refer to the placement of poles and zeros in the z-domain helping us achieve convolution, if FIR and ??? in IIR.

Chebyshev and Butterworth define certain constraints under which a filter can be designed and is indicative of the ripple effect in the pass-band and stop band.

My question is - how are these two related. Can we design an IIR using the Chebyshev constraint. In real-time, is this constraint used using the Impulse Invariance Method or Bilinear Transform for IIR design?

Any answer to this would really help. Thanks

Chebychev and Butterworth filters can be transformed with the bilinear transform into IIR filters.  The amplitude and phase responses will be warped according to the bilinear transform.  This warp mainly makes for a higher cutoff rate than the corresponding analog filter (handwavingly, it's because the antialiasing filtering and sampling already does part of the work) which tends to be ok for most applications.

In contrast, using a bilinear transform on a Bessel filter is quite pointless since its salient design feature is a near constant phase delay in the pass band, and the bilinear transform warp messes this up.  So there is a separate digital equivalent to Bessel filters that isn't related via simple bilinear transform.

FIR filters, in contrast, don't have straightforward connections to analog filter designs.  Their design procedures are quite separate and only suitable for digital filters.


Spectrum and Spectrogram explanation I'm reading this material but I don't understand this example:

"An example of a 94% overlap transform processing could be demonstrated by a 279 sample frame with 262 samples overlapped from each frame into the next. At 12.5 MHz sample rate, each spectrum frame would occupy 23.76 µs. Each spectrum would start 1.36µs after the previous one and could show a 10 MHz span."
Could anybody please help to explain what is 10 MHz span and why why each spectrum could show a 10 MHz span?
The same problem with the following picture (I don't understand the example):

Here's an attempt to answer your first question:

  "An example of a 94% overlap transform processing could be demonstrated by a 279 sample frame with 262 samples overlapped from each frame into the next.

A spectrogram is really just a sequence of spectra.  In this case, each individual spectrum is taken of 279 time samples (which I will call a frame).  Each subsequent frame consists of the last 262 samples of the previous frame, plus the next 17 time samples.

That means percentage of common samples between subsequent frames is:

$$
\frac{262}{279} = 93.9068 \%
$$

  At 12.5 MHz sample rate, each spectrum frame would occupy 23.76 µs.

We have 279 samples in each frame.  The distance between time samples is

$$
\Delta t = \frac{1}{12.5 \times 10^6} = 80 {\rm ns}
$$

Add 279 of those together and we get

$$
279 \times 80 {\rm ns} = 22.32 {\rm \mu s}
$$

It looks like someone can't multiply:  the number 23.76 µs is actually $297 \times 80 {\rm ns}$.

  Each spectrum would start 1.36µs after the previous one...

Each frame starts 17 samples after the last.  This is $17 \times 80 {\rm ns} = 1.36 {\rm \mu s}$

  ... and could show a 10 MHz span." Could anybody please help to explain what is 10 MHz span and why why each spectrum could show a 10 MHz span?

It's not clear what they are meaning by this.  If you sample at 12.5MHz, then the total possible frequency span is 12.5MHz (if you include positive and negative frequencies).  There may be some allowance for guard bands that reduces the usable bandwidth from the ideal.


Filtered signal does not follow my cut off frequency I filtered my signal with low pass filter butterworth forth order with cutoff frequency of 50 hz. However when I compare the filtered signal with original one in frequency domain, it seems that the cut off frequency started from 25 hz instead of 50.

this is my matlab code and my data from here

%sampling frequency is 146

[b,aa] = butter(10, 50/146,'low');
data_filtered = filtfilt (b,aa,data);

did i make some mistake ?

did i make some mistake ?

Yes. The cutoff frequency is specified relative to the Nyquist Frequency and not the sample rate

From the documentation

  Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample
  rate.

Also: you design a 10th order butterworth and filtfilt() doubles that again to 20. That will give you significant domain "ringing".


How to show that the signal $x_n = A\cos(\omega n)$ can be fully predicted by a system with two weights $w_1,w_2$ I am trying to solve the following exercise:

  Show that the signal $x_n = A\cos(\omega n)$ can be fully predicted by a system with two weights $w_1,w_2$ (i.e. $x_n = w_1 x_{n-1} + w_2 x_{n-2}$). Find $w_1,w_2$.

Some ideas came up but, even if they are right, I can't seem to order them in the right direction:

I can see that $x'' = -\omega^2x$. Does it connected somehow to what is asked?
I guess it has something to do with stationary signals. If yes, how can I prove that this signal is stationary? (Also, somehow I can't find in the internet a clear definition of "stationary signal". Just vague ideas like "not dependant on time". What does this mean mathematically for a signal to be "not dependant on time"?)
Maybe it is concerned with some trigonometry of this form:
\begin{align} x_{n+2} &= A\cos((n+2)\omega) \\ &= A\cos(n\omega + 2\omega) \\ &= A\cos(n\omega)\cos(2\omega) - A\sin(n\omega)\sin(2\omega) \\ &= B\cos(n\omega) + C\sin(n\omega)\end{align}
I could go on with this development but I am realy not sure what am I looking for..

Any help?
Thanks!

HINT:

Use

$$2\cos(x)\cos(y)=\cos(x+y)+\cos(x-y)$$

with $x=(n-1)\omega$ and $y=\omega$.


Difference equation with variable coefficients in MATLAB Given the following difference equation

$$y[n]-\left(\dfrac{n}{n+1}\right)y[n-1] = n\cdot x[n]$$

How can we use MATLAB to solve it?

I know if the coefficients are constant we can simply use filter(b, a,x), but not sure how to do it with variable coefficients.

When $n=0$ you may need to know $y[-1]$. The current output $y(n)$ depends on the current input $x(n)$ and previous output $y(n-1)$ scaled.


Most accurate way to find RMS of signals fundamental I would like to know which algorithmn is the most accurate one in order to calculate the fundamental RMS of two signals (f. e. voltage and current). Is it the standard FFT approach with a rectangular window or may there be some more precise algorithmn?

In this example, the resulting RMS signals may be used further for calculating the signal of the apparent power. The following Python code may help to illustrate the question:

import numpy as np

#sampling frequency in Hz
Fs = 16000

#fundamental frequency in Hz
F1 = 50

#number of samples
samples = 1000000

#generate voltage and current signals
t = np.arange(samples / Fs)
u = np.sin(F1 * 2.0 * np.pi * t) + 0.5 * np.sin(80.0 * 2.0 * np.pi * t) + 0.1 * np.sin(10.0 * 2.0 * np.pi * t) #signal may have additional harmonic and sub-harmonic content
i = np.sin(F1 * 2.0 * np.pi * t) + 0.2 * np.sin(120.0 * 2.0 * np.pi * t) + 0.3 * np.sin(30.0 * 2.0 * np.pi * t) #signal may have additional harmonic and sub-harmonic content

#apply some algorithm to find accurate RMS of the signals fundamental
#t_rms_50Hz = ...
#U_rms_50Hz = ...
#I_rms_50Hz = ...

#get apparent power S
S_50Hz = U_rms_50Hz * I_rms_50Hz

This is the start of an answer containing my initial thoughts that may be helpful (I wrote quickly due to limited time so may have errors, I will review in more detail later but can at least offer some direction):

One approach is to use this resonator that was developed in this post to isolate your frequency of interest: FFT analysis for Vibration Signal. The issue you will have is allowing for enough time for the filter to settle to it's final amplitude on the output (which is based on the bandwidth of the filter, translated to time constants in time, and depending on your accuracy needed is how many taus you will need to wait until you have signal that you can take the rms value of. I have not confirmed this but it also seems that you could be clever and use the known information of settling time so that you can utilize the initial settling signal's information as part of your answer (weighting the rms calculation by the result in an "optimal combining" fashion.

You could also do FFT's with windowing in which case refer to this paper by fred harris regarding the effect the window will have on your rms calculation: fred harris On the Use of Windows

Note with an FFT approach that each bin in the FFT is a filter that is responsive to the signal anywhere within its filtering bandwidth. Without any further windowing, it is a "rectangular" window which has the highest frequency selectivity but lowest dynamic range, as evidenced in the plot below of a simple 4 point DFT. Any signal power that "leaks" into the bandwidth of an FFT bin will become part of your RMS calculation. The filter responses as shown below approach Sinc functions as the number of points increase, with the null to null spacing of the main lobe for each filter as 1/N where N is the number of points in the FFT. When you use a window, the sidelobes significantly decrease, but at the expense of the main lobe.

Given you have spurs that you are concerned with that are only 6 Hz away from a signal that is 16KHz, if you assume you used a window that increased your main lobe 3x (for example) to get your desired sidelobe attuenuation, meaning was 3/N null to null, Then you would need a minmumum FFT length on the order of $3\times 16000/6=8000$ which is not so bad.  (you can take the FT of various windows to see for yourself what each would do- my favorite is the Kaiser window because you can adjust selectivity and bandwidth with the additional $\beta$ factor).

Note these particular items that will effect you if using a window: coherent gain, noise gain, scalloping loss.

Coherent Gain: The standard deviation of a single tone will be modified by the coherent gain of the window. The coherent gain is found by summing all the samples in the window itself:

$$G_c=\sum_0^{N-1}w_n $$

For instance, the coherent gain of a rectangular window is N, so the net gain of a window compared to the rectangular window is $G_c/N$, which if you try yourself with common windows, you will matches what is tabulated in table 1 in fred harris' paper.

The noise gain however is how the standard deviation of the noise components in your signal will scale, where noise is anything that is uncorrelated from sample to sample. The noise gain due to a window increases as the sum of the squares:

$$G_n=\sqrt{\sum_0^{N-1}w_n^2}$$

So note the processing gain that results, for the rectangular window the signal will increase by $20log_{10}(N)$ while the noise will increase by $10log_{10}(N)$ for a net gain of $10log_{10}(N)$ in SNR. You can use this and the measurement of the variance of your entire signal for a more accurate estimate of the signal component since the variance measurement is S+N, and your post FFT measurement after scaling will be S+N/K where K is the resulting processing gain for your window used.

The final item to mention is "scalloping loss" which occurs if your signal of interest is not in the center of the bin (observe the immediate roll-off of a sinc function). You can either account for this mathematically in your coherent gain if you have specific knowledge of your signal, or you can minimize the effect of scalloping loss by using overlap-add techniques (detailed in fred harris' paper).

Having written this out, between the FFT appraoch and my first suggestion of using the tuned 2nd order resonator, I favor the FFT approach as being simplest and most straightforward. You also have significantly more samples than the minimum needed to reject your interference signals, which are sufficiently low so you should be able to extract a good estimate of the standard deviation of your signal. Note the overall SNR is also of interest as you can use that to establish the confidence interval for your result. For this see Estimating confidence intervals for the mean value


Where can we find noises like Gaussian, Rayleigh, Gamma, exponential noises in real world images? I know that while taking pictures from a camera, noise may appear in the image.
It is often described as  Gaussian, but what about the other noises, when do they occur in real world?

Please help!

Gaussian Noise: Standard model of amplifier noise.

Gamma Noise: Frequent in Radar/LiDAR.

Rayleigh Noise: One example where the Rayleigh distribution naturally arises is when wind velocity is analyzed into its orthogonal 2-dimensional vector components.

Exponential Noise: Channel Based Communication(see AEN).


Fourier coefficients of odd and even part of a signal I have this signal  and I have to find the Fourier coefficients of the odd and even part. First I found that $$ x_p(t) = \frac{1}{2} ( x(t) + x(-t) ) $$ and I made the graphic of this part and I obtained a straight line in 1/2. I made the same thing for $$ x_d (t) = \frac{1}{2} ( x(t) - x(-t) ) $$ and making the graphic I obtained the signal $$\mathrm{rect}(\frac{ t- n \frac{T_0}{4} }{ \frac{T_0}{2} }) - \frac{1}{2} $$.

Now $$X_{pk} = \frac{1}{T_0} \int_{0}^{T_0} \frac{1}{2} e^{-j 2 \pi k f_0 t } dt $$  and I obtained that  if $$ k=1 $$ $$X_{pk} = 0 $$ and if $$  k=0 $$ $$X_{pk} = \frac{1}{2} $$. This result is the same of my book. For the odd part I should do $$X_{dk} = \int_{0}^{T_0}( \mathrm{rect}(\frac{ t- n \frac{T_0}{4} }{ \frac{T_0}{2} }) - \frac{1}{2})e^{-j 2 \pi k f_0 t } $$ and I obtain $0$ if $k=1$ and $1$ if $k=0$ but the result should be $0$ if $k=0$ and $$ \frac{-j}{k \pi } $$ if $k = 1$.

OK, so your signal is described in one period, $T_0$, as $$x_{T_0}(t) = \left\{\begin{array}{ll}
0, & -T_0/2 < t \leq 0 \\
1, & 0 < t \leq T_0/2
\end{array}\right.$$
Let's keep it that way (without any other shortcuts like $\prod$ or $\mathrm{rect(\cdot)}$ etc).
One way to find the Fourier coefficients of the even and the odd part of a signal is to actually find the even and odd part of the signal, and do a Fourier Series analysis on each.

But, another way is to check Fourier Series properties - again :-) .

If $x(t)$ is real - which it is - with Fourier coefficients $X_k$, then the odd part of $x(t)$, let it be $x_{odd}(t)$, has Fourier coefficients $X^{odd}_k = j\Im\{X_k\}$, where $\Im\{X_k\}$ is the imaginary part of the Fourier coefficients of the original signal $x(t)$.

Similarly, the even part of the signal $x(t)$, let it be $x_{ev}(t)$, has Fourier coefficients $X^{ev}_k = \Re\{X_k\}$, where $\Re\{X_k\}$ is the real part of the original Fourier coefficients $X_k$.

In short: $$x(t) \longleftrightarrow X_k$$ $$x_{odd}(t) \longleftrightarrow j\Im\{X_k\}$$ $$x_{even}(t) \longleftrightarrow \Re\{X_k\}$$

So, you have two options:

Find $x_{odd}(t)$, $x_{ev}(t)$ and find their Fourier coefficients analytically (takes a lot of time and paper)
Find the Fourier coefficients, $X_k$, of $x(t)$ and then take their real and imaginary part (faster and easier)

I would go for the second option, since $x(t)$ is very simple and its Fourier coefficients are given by $$X_k = \frac{1}{T_0}\int_0^{T_0}x_{T_0}(t)e^{-j2\pi kf_0t}dt = \frac{1}{T_0}\int_0^{T_0/2}1\cdot e^{-j2\pi kf_0t}dt$$

P.S: In case you wish to go with the first option, the even part of the signal can be computed by first finding $$x_{T_0}(-t) = \left\{\begin{array}{ll}
1, & -T_0/2 \leq t \leq 0 \\
0, & -T_0 < t < -T_0/2
\end{array}\right.$$
and then $$x_{even, T_0}(t) = \frac{1}{2}(x_{T_0}(t) + x_{T_0}(-t)) = \frac{1}{2}$$

Similarly, the odd part can be obtained as $$x_{odd,T_0}(t) = \frac{1}{2}(x_{T_0}(t)-x_{T_0}(-t)) = \left\{\begin{array}{ll}
-\frac{1}{2}, & -T_0/2 <t < 0 \\
\frac{1}{2}, & 0 \leq t \leq T_0/2
\end{array}\right.$$

So the even part is just the mean value of the signal, $X_0$. The odd part of the signal has Fourier coefficients $$X_k = -\frac{1}{T_0}\int_{-T_0/2}^0 \frac{1}{2}e^{-j2\pi kf_0t}dt + \frac{1}{T_0}\int_0^{T_0/2}\frac{1}{2}e^{-j2\pi kf_0t}dt$$

I think you can continue from here. Let me know if you make it till the end. :)


Filter to eliminate a 60 Hz frequency range I have a analog signal that is converted into a discrete-time signal with an ideal A/D converter with a sampling frequency $fs$.  The bandwidth of the signal of interest is $1 kHz$. The resulting signal $x[n] = x_a(nT_s)$ is then processed with a discrete-time system that is described by the difference equation:
$$y[n] = x[n]+ax[n-1]+bx[n-2]$$
The filtered signal, $y[n]$, is then converted back into an analog signal using an ideal D/A converter. I need to determine values for $fs$, $a$ and $b$ in order to remove a $60Hz$ interference that has a signal in the form: $i_a(t) = Asin(120 \pi t)$
By Nyquist, I know that the value of $fs$ needs to be greater than $2000Hz$, however, I don't know how to determine $a$ and $b$. Any thoughts?

The OP is trying to implement a notch filter with a 3 tap FIR filter (as restricted by the problem he is trying to solve). It must be noted that this would result in a very poor filter implementation: It will reject 60 Hz perfectly but will have unavoidable attenuation over a very broad section of the desired passband range. An actual practical implementation would be done with an IIR approach in order to provide a tight notch, as demonstrated by this link: Transfer function of second order notch filter
As a hint to determine the filter within this FIR restriction, as an exercise only and of no real practical value proceed as follows:
Take the Z transform of the difference equation.
Determine the transfer function form as a polynomial ratio $H(z) = Y(z)/X(z)$
From this know that the frequency response is determined by restricting $z$ to the unit circle ($z= e^{j\omega}$) with $\omega$ extending uniquely over the normalized radian frequency from $-\pi$ to $+\pi$ which corresponds to $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate.
Given this is an FIR structure, the transfer function will result in simply a 2nd order polynomial in the numerator (whos roots are the "zeros" of the filter). Place the zeros of this polynomial such that they corresponds to the 60 Hz locations on the unit circle (+/- 60 Hz for a real filter, mapped to the angular frequency as described above), which would provide complete attenuation at f = 60 Hz.


Wiener Filter in Frequency Domain: What it does to a specific Frequency? As I understand Wiener filter in time domain tries to estimate a signal as close as possible to its (original) non-degraded signal using the degraded signal by white noise.

$$H(\omega)=\frac{\Phi_{ss}(\omega)}{\Phi_{ss}(\omega)+\Phi_{nn}(\omega)}$$

How does this function work in the frequency domain? Say for an example at 30 Hz there is noise and how the transfer function looks like?

You already know the formula, so that's how the frequency response looks like. In frequency regions where the signal is stronger than the noise, the frequency response is close to $1$, so the portion of the input signal in that frequency range passes the filter with almost no attenuation, whereas in frequency regions where the noise is stronger than the signal, the filter's frequency response becomes small, and the corresponding part of the input signal is attenuated.

If you define a frequency dependent SNR

$$\textrm{SNR}(\omega)=\frac{\Phi_{ss}(\omega)}{\Phi_{nn}(\omega)}\tag{1}$$

you can rewrite the frequency response of the Wiener filter as

$$H(\omega)=\frac{\textrm{SNR}(\omega)}{1+\textrm{SNR}(\omega)}\tag{2}$$

Eq. $(2)$ clearly shows that if $\textrm{SNR}(\omega)\gg 1$ you have $H(\omega)\approx 1$ (i.e., no attenuation), whereas if $\textrm{SNR}(\omega)\ll 1$ you get $H(\omega)\approx \textrm{SNR}(\omega)\ll 1$ (i.e., strong attenuation).


Best way to inverse Hilbert on FPGA? I have a real signal that is being sampled via ADC and passed to a Xilinx V6 FPGA.  I then do a Hilbert transform on it using the FIR compiler core and get the complex signal (I and Q) to tweak the signal some.  Now I am struggling with how to do the inverse Hilbert to get back to a real signal again.

I understand in theory that H(H(u)) = -u.  But the FIR compiler simply does not allow for a complex input and a real output.  Is there a trick to get around this?

Just take the real part of the complex signal.

(And if you do not believe it: click here.)


What would be the outcome of the convolution of two images? Lena

Mask

I have written a C# code which produces the following output.

Is my output correct?
Am I padding in the correct way?

You can use MATLAB to get the 2D convolution of these images. Before you do so, you'll want to convert image 2 to grayscale and convert both to double.

The result has a very large range, so I normalized the second image.

lena = im2double(imread('lena.jpg'));
color = im2double(rgb2gray(imread('color.jpg')));
out = conv2(lena, color./sum(sum(color)));

If you happen to be interested in blending images together, you might want to check out this paper.


Good book or reference to learn Kalman Filter I am totally new to the Kalman filter. I've had some basic courses on conditional probability and linear algebra. Can someone suggest a good book or any resource on the web which can help me can understand Kalman Filter operation?

Most websites start directly with the formula and what they mean, but I am more interested in its derivation, or if not detail derivation then at least the physical significance of each operation and parameter.

Many years ago I wrote this tutorial on the Kalman filter.  It derives the filter using both the conventional matrix approach as well as showing it's statistical assumptions as an 'optimal' least squares filter.


What transform can I use to "zoom" in a specific portion of a signal's spectrum? If I recall correctly, there is a variation of the DFT that can be used to analyze a specific band of the spectrum of a signal. How is it called?

You're probably referring to the zoom FFT. It's essentially a technique that allows for complexity reduction in the case where you have a small portion of a larger band that you'd like to analyze at high spectral resolution. It prevents you from having to calculate the high-resolution frequency content in the bands that you don't care about. Roughly, the algorithm can be summarized as follows:

Apply a bandpass filter around the region of interest, eliminating the components outside the band that you care about.
Decimate the signal by a factor $D$, such that the resulting sample rate still meets the Nyquist criterion for the filter's passband width. Depending upon where the band's center frequency was, this process might also involve frequency-translating the signal to baseband.
Perform a DFT on the signal. In order to get the same frequency resolution at the output, the "zoomed" transform only requires a transform length that is $\frac{1}{D}$-th of what you would need to use if you used the original, unfiltered signal for your analysis.


Extracting motion vectors from H.264 bitstream I'm looking for an open-source tool/code or some guidance to extract the motion vectors (MVs) of a H.264 encoded bit sequence. I'm already aware that motion vectors can be visualized using ffmpeg with the following command:

ffplay -flags2 +export_mvs input.mp4 -vf codecview=mv=pf+bf+bb

However, I want to produce a log file where the MVs of P and B frames are listed frame by frame. I checked out the structure of MVs from libavutil/motion_vector.h, but I couldn't find an example which shows how they are extracted and laid over the original sequence by ffplay. I thought that if I can find that out, I could possibly re-arrange the code to extract the MVs to a text file.

I would appreciate any example codes or hints.

I think the following code can help you:
https://github.com/vadimkantorov/mpegflow
The main idea is that you need to decode frames first and query motion frames by calling av_frame_get_side_data function.


Implementation of Hilbert Spectrum in MATLAB I am trying to implement this equation  on Matlab,  where i have wj values in the form of an array, cumtrapz command of MATLAB implements this, but its confusing which way is exact. such as cumtrapz(t,wj) or simply cumtrapz(wj) or manually doing it? Because results of all these though related, are different.

If $\omega$ is a vector with all the discrete elements $\omega_{j}$ over the descrete time steps $dt$ the correct formula would be dt*cumtrapz(omega); which is the same as cumtrapz(dt,omega) or alternatively if $dt$ isn't constant but $t$ is a vector with evaluation positions Omega=cumtrapz(t,omega,2). The "2" is needed to indicate that the 2nd vector should be integrated.


Properly interpreting FFT output I understand that the DFT (Discrete Fourier Transform) process produces the same number of outputs as there are inputs, and this is clear from the basic DFT expression

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-i2\pi kn/N} \tag{1}$$

However, when computing the FFT (Fast Fourier Transform), the final output is the sum of two N/2 transforms (the even and odd DFTs).

$$X[k]=X_{even}[k]+W_{N}^{k}\cdot X_{odd}[k] \tag{2}$$

Suppose $N=16$. If I try to evaluate Equation (2) for $N=12$, it won't work, since $X_{even}[12]$ doesn't exist (since $X_{even}$ only has $N/2=8$ inputs in the first place).

This is normally solved, as in the case of Fat32's solution to a previous question of mine, by concatenating the outputs of each $N/2$ DFT with themselves to produce an output that is N long, which would result in an $X_{even}[12]$ which does exist. Alternatively, one could add a condition that if $k>N/2$, then substrat $N/2$ from each $X_{even}$ and $X_{odd}$ index.

I was wondering what the mathematical foundations of each of these work around are, since they aren't explicitly included in expression 2? Is it something that is basically common sense? Or does it really mean that only half of the outputs are actually valid?

I also have a hard time interpreting the output of the FFT since the second half of the outputs are mirrored and actually represent negative frequencies. So doesn't that mean that in order to really present an accurate output I must chop that second half off and move it to the negative frequency domain?

There seem to be several ways to present the output (some leave the negative frequencies mirrored and in the positive domain, and some do order them properly). Any advice about the best way to interpret and present the outputs of the FFT process is greatly appreciated!

The upper half of the DFT being the conjugate mirror of the lower half is only true for real valued signals.

For an alternative understanding of what the DFT means you can read my blog article DFT Graphical Interpretation: Centroids of Weighted Roots of Unity.  I also recommend your read my first article "The Exponential Nature of the Complex Unit Circle" for an understanding of the foundation of the DFT, in particular, the significance of the Roots of Unity.

Hope this helps.

Ced


Is the calculation of variance correct from an estimator -- confusion regarding complex number where I am going wrong? I have an expression for the varaince of measurement noise obtained from an estimator. The measurement noise is additive white gaussian having values in complex domain. I have squared the term $(.)^2$, should I be squaring or just implement without squaring? But if I dont square the term, then I am getting a complex valued variance which is incorrect as variance is always real valued.
Kindly suggest what to do.

You must put:

$$\sigma^2=\frac{1}{N-1}\sum_{i=1}^N(y_i-\hat{y}_i)\overline{(y_i-\hat{y}_i)}$$

that is, replacing your code by:

for i = 1: N
   dy(i,1) = (y(i)-yhat(i))*conj(y(i)-yhat(i));
end
sigma2 = 1/(N-1)*sum(dy);

The variance is the second moment of a probability distribution, which can be defined for any metric space in the sense of the definition over the random variable (integral of a density), or in any vector space, in the sense of the estimator of the random variable (sum of n samples).

For $X$ a random variable  $\in \mathbb{R}^n$, with probability density function $p(x):I_X \subset \mathbb{R}^n \rightarrow \mathbb{R}$, the second moment $\in \mathbb{R}$ is:

$$\sigma^2=\mathbf{E}\{(X-\mu)^2\}=\int_{I_X} (x-\mu)^2 p(x)dx$$

which is the variance for $n=1$ and the covariance matrix in general.

For other fields $\mathbb{F}$, complex and even real numbers, the metric $d(x,y): \mathbb{F}x\mathbb{F}\rightarrow \mathbb{R}$ can redefine the moment definition

$$\sigma^2_{\mathbb{F}}=\mathbf{E}_{\mathbb{F}}\{(X-\mu)^2\}=\int_{I_X} d(x,\mu)^2 p(x)dx$$

For example, for the complex field, the metric function can be the standard distance from the hermitian inner product, and thus we have the variance definition:

$$d(x,y)=\sqrt{(x-y)\overline{(x-y)}}$$
$$\sigma^2_{\mathbb{C}}=\mathbf{E}_{\mathbb{C}}\{(X-\mu)^2\}=\int_{I_X} (x-y)\overline{(x-y)} p(x)dx$$

But the standard product, without the conjugated term, though can be a valuable figure (named as the "complementary variance") is not a second moment:

$${\sigma^2}_{\mathbb{C}comp}=\int_{I_X} (x-y)^2 p(x)dx$$

Of course, you always could be more "reasonable" and prefer to drop the complex analysis and move onto a $\mathbb{C}^2$ cartesian plane, and start taking covariances instead of a single number which dont dive you the full distribution information of your stochastics.


How does MATLAB handle IIR filters? MATLAB has a butter, which constructs a Butterworth filter given an order and relative cutoff frequency. The filter created can be used to filter any finite signal.

How does MATLAB do this if the filter has an infinite impulse response? I suppose it must window the signal -- is this done by a simple rectangular window so that all known values are used?

Also, how is the continuous filter discretized? Bilinear transform, impulse response matching, etc?

I don't think the MATLAB help text explains this:

  Y = FILTER(B,A,X) filters the data in vector X with the
     filter described by vectors A and B to create the filtered
     data Y.  The filter is a "Direct Form II Transposed"
     implementation of the standard difference equation:

a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
                      - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

This isn't really a MATLAB-specific issue; I see a couple more general questions:

How do you implement a digital IIR filter?

You can apply any general digital filter by convolving its impulse response with the signal that you want to filter. That looks like:

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

This works great for FIR filters, but you run into complications with IIR systems, because their impulse responses are infinitely long ($N \rightarrow \infty$ in the above sum). That makes the filter output difficult to calculate using the straightforward convolution sum above. To actually implement an IIR filter, we need a way to express the filter output in a way that is more computationally tractable.

The attribute of an IIR filter that makes its impulse response infinitely long is the system's recursive nature; there is feedback from the output back to the filter input. This means that the output of an IIR filter at any given instant of time can depend upon two separate information sources:

The input signal to the filter, past and present, and
Past values of the filter output signal.

This leads to the difference equation representation of the IIR system:

$$
a_0 y[n] = \sum_{k=1}^{M-1} a_k y[n-k] + \sum_{k=0}^{N-1} b_k x[n-k]
$$

This representation illustrates that we can implement an IIR filter by calculating a weighted sum of $M$ past outputs of the filter and $N$ past inputs of the filter (in practice, $N$ and $M$ could and often are the same; the order of the filter when defined as above is $max(M,N) - 1$). This is a closed-form expression that captures the behavior of the filter fully, suitable for automated implementation, and is the equation referenced in the MATLAB help excerpt that you showed in your question. Note that $a_0$ is almost exclusively assumed to be equal to $1$; I only included it because it is shown in the quoted help text, referenced as a(1).

So to specifically answer your question, as the text suggests, MATLAB implements IIR systems using the above equation. When you design an IIR filter with MATLAB or otherwise, you get two sets of coefficients ($a_k$ and $b_k$ in the difference equation) that define how to weight past filter inputs and outputs in the filter implementation. There are some other subtleties that are sometimes relevant; for instance, the text refers to using a "direct-form II transposed" realization of the filter. You then apply the filter by calculating the difference equation above for each desired output sample.

There are in fact a number of topologies that can be used to implement digital filters, each with their own tradeoffs. DF2T is often used because of its more efficient structure (it contains a minimum number of delay elements) and its improved robustness to roundoff error (since the feedforward coefficients $b_k$ are implemented first).

How does MATLAB design Butterworth filters?

The documentation for the butter function does not specify how it generates the discrete approximation to the analog Butterworth filter prototype. However, by peeking into butter.m, you find:

% step 5: Use Bilinear transformation to find discrete equivalent:
if ~analog,
    [a,b,c,d] = bilinear(a,b,c,d,fs);
end

So it seems that it uses the bilinear transform to map the Butterworth filter prototype to a digital filter realization.


How to generate sinusoidal images that shift over time I try to generate an image using this code.

w=10*pi/256;
E=zeros(256,256,100);
for t=1:100
for y=1:size(E,1)
    for x=1:size(E,2)
        E(x,y,t)=255*(sin(w)*sin(w*y*2)+1)/2)+1;
    end
end
imshow(E(:,:,t),[])
title(t);
drawnow
pause(0.1)
end

However I'm confused how to generate the image that shift over time, like for time t=2, the stripes is shifted by 0.0001 in the x direction.

Any idea how to fix this? Thank you.

Your code will produce a constant pattern along $x$ direction, so shifting this constant pattern along $x$ direction won't change anything.

To shift along $y$ direction you could replace $sin(2wy)$ by $sin(2wy + ct)$, and you have to choose $c$ depending on the amount of shift that you want.

Also you could create an initial image (for $t=1$) with an exactly integer number of periods, if at each time step your image shifted by integer number of pixels use circshift function to shift it by the amount that you want. If it doesn't shift by an integer number of pixels, you  have to take the FFT of the image and change the phase linearly. in this case you have to change the phase by $c\times W_x$ where $c$ is the ammount of shift in the $x$ direction and the $W_x$ is the spatial frequency in the $x$ direction.


Can someone explain me how to generatee a doppler spectrum? I have a signal generated in matlab, with Rayleigh Fading, how can I add a given type of Doppler spectrum ( Gaussian for instance) to that signal ?

Any help will be really appreciated

When we talk about Doppler, we typically refer to time-varying systems.

A popular model in system theory is to assume systems to be LTI (linear time-invariant). This means that beside being linear, their input-output transformation does not change in time, it does not matter when we apply a certain input signal $x(t)$, we'll always get the same $y(t)$. Such systems are completely characterized by an impulse response $h(\tau)$ which specifies exactly this transformation: $y(t) = h(t) * x(t)$, where $*$ denotes convolution.

In practice, this assumption is often not true. Wireless communications are a good example. While the channel is typically (approximately) linear, its characteristics change in time, e.g., due to movement of transmitter, receiver or interacting objects such as cars, trees, people, etc. For this reason, we need to extend the impulse response by another time variable that tracks the change of $h(\tau)$ in time. We often simply write $h(\tau,t)$. It may seem weird to have two time variables but there is solid theory behind this. As an engineering perspective, we can think of $\tau$ as the "fast time" (at which signals travel, orders of $\mu$s) and $t$ as the "slow time" (at which the channel changes, order of ms) and treat them more or less as independent.

The rate at which the channel changes is called Doppler frequency, say $\alpha$. Mathematically we can think of it as the Fourier transform of the time variable $t$. Given that we have two time variables and we can Fourier both of them, this leads to four possible representations of the input-output relation: $h(\tau,t)$ (time-varying impulse response), $H(f,t)$ (time-varying transfer function), $h(\tau,\alpha)$ (Delay-Doppler spread function), and $H(f,\alpha)$ (name it as you like). They are also referred to as the (first set of) Bello functions to drop a keyword for further reading.

Since this is the title of your question: the Doppler spectrum measures how much power is associated to a certain Doppler frequency, it is the power spectral density of $h(\tau,t)$ over the time variable $t$.

Regarding the implementation of a certain Doppler spectrum into your simulations: in general what you would do is to generate the channels for different time points independently and then add a filter over time that introduces the proper amount of correlation. If you're using Matlab, have a look at [MATLAB Fading channel page], it explains it quite well.