Butterworth filter approximation: derivation and output poles I am having trouble understanding the exact derivation of the butterworth filter and how it results in the output of the poles. I have researched multiple lecture series and textbooks and this is my understanding so far:

An idealised low-pass filter, i.e. a brick wall, can be realised with the equation below as as $n$ tends to infinity?
$$
\lvert H(j\omega)\rvert = \frac{1}{\sqrt{1+\omega^{2n}}}
$$

This is the amplitude response - i.e. the magnitude of the frequency response?
Now by squaring this amplitude response and setting $s = j\omega$ you get the following:
$$
\lvert H(s)\rvert^2 = \frac{1}{1+\left(\frac sj\right)^{2n}}
$$
What does squaring the magnitude of the frequency response get you?
I understand that $\lvert X(f)\rvert^2 = X(f)\cdot X^*(f)$?
The result of these two terms(/transfer functions?) $H(s)$ and $H(-s)$. Subsequently you are only interested in the $H(s)$, why? I understand it's poles lie on the left of $s$-plane and therefore the system is stable? but what about $H(-s)$?
You then obtain the poles from the resulting derivations and get a resulting transfer function. This transfer function and its related co-efficients can be implemented physically with the appropriate R and C values?

Any help would be greatly appreciated!

Regards

Note that the Butterworth filter is only one of many possible approximations of an ideal low pass filter. It is a Taylor series approximation (at $\omega=0$) of the desired response, and, consequently, the magnitude response of the Butterworth filter is maximally flat at $\omega=0$.

Approximating the squared magnitude response (instead of the complex frequency response, or the magnitude response) has the advantage that the function to be approximated is an analytic real-valued function of a real-valued variable, which simplifies the approximation problem.

We have

$$|H(j\omega)|^2=H(j\omega)H^*(j\omega)=H(j\omega)H(-j\omega)=H(s)H(-s)\Big |_{s=j\omega}\tag{1}$$

because $H^*(j\omega)=H(-j\omega)$ for real-valued systems, i.e., for systems with a real-valued impulse response, or, equivalently, with real-valued transfer function coefficients.

Since for a Butterworth low pass filter

$$|H(j\omega)|^2=\frac{1}{1+\omega^{2N}}\tag{2}$$

we have

$$H(s)H(-s)=\frac{1}{1+(s/j)^{2N}}=\frac{1}{1+(-s)^{2N}}\tag{3}$$

From the squared magnitude function we need to derive the transfer function $H(s)$, which is the function that we actually want to implement. From $(3)$ we know that the poles of $H(s)H(-s)$ (i.e., the zeros of the denominator of $(3)$) lie one a circle in the complex plane. Since we can only implement a causal and stable transfer function, we choose the $N$ poles in the left half-plane for $H(s)$, and consequently, the mirrored poles in the right half-plane are the poles of $H(-s)$.

Now we know $H(s)$ (i.e., its numerator and denominator coefficients), and, after having chosen a cut-off frequency and a filter structure, we can compute the values of the passive elements necessary to implement the given transfer function.


Effectively extracting “real” data from a noisy dataset This question was originally posed to Math SE. It was suggested that DSP SE would be more suitable.

Background & Motivation:
I have three lists of timestamps (UNIX timestamps, plus a subsecond part (e.g. 1451606400.9475839201)). They are not necessarily of the same size, however they are ordered. Each such list corresponds to a real-life instrument, of which there are three. The instruments in question record a time stamp each time they observe an "event". We may assume that their clocks are synchronized.
The issue is that the instruments are so sensitive, they record timestamps very frequently (on the order of 10 Hz), and only a small portion correspond to actual events. It is difficult to put a number on precisely how many timestamps are real events (perhaps a handful per year).
Further, the instruments are close in proximity, so if one instrument "sees" and event, the other two almost certainly do as well. The clocks are precise enough to record the difference in observation time, even over miniscule differences.
The data is for one year, the same year for each of the three lists. We may assume that the "random" timestamps are uniformly distributed. There are, however, gaps in the data (e.g. maybe all three instruments were off for the months of March, April, May). The gaps will be the same for all three lists.

Goal:
Using only the timestamps, I want to attempt to find those which are "likely" to correspond to a "real" event, such that further analysis can be conducted. The "events" in question are light signals, so I can restrict my search to only those for which the difference in observation for some triplet of timestamps is less than the transit time between two observers.
My first inclination was to, from the three lists, produce a list of triplets, one timestamp contributed by each list, such that $max(A,B,C) - min(A,B,C)$ was minimized. Unfortunately, this found very few "triplets" which a) fell in the restriction mentioned previously and b) corresponded to "real" events. I mentioned previously that there are very few events, however I expect there to be far more than what this found.
I then tried doing the above, but minimizing the $\chi^2$ error, which I defined for some triplet $A,B,C$ (one from each list/instrument) as $(A-B)^2 + (A-C)^2 + (B-C)^2$. This found even fewer "triplets", and no real events.

Problem:
What techniques can be used to extract "coincidences" (the "real events") from a set of more-or-less uniformly distributed data, where there is far more random than "real" data? Here, we assume that I have access only to the timestamps.

Judging from your description alone, I'd say an algorithm to find appropriate 3-clusters of timestamps would be something like:

compare the three first elements of your sorted lists. Pick the earliest of the three.
Remove the element from its list.
determine difference of that element to the other two first

If both difference are below a "simultaneity" threshold, register a new "common event".

return to 1

Often, you'd also remove the two "matches" from the other two list in 3., but that reduces your ability to detect ongoing physical excitation of your sensor as one by post-processing your "common events" list.
Seeing you're only dealing with about a billion timestamps in total, this simplified variation of your distance-limited cluster detector should allow for quick experimentation¹.

But:
Reading through your question, it seems likely that one physical phenomenon can lead to multiple adjacent timestamps in a single list, so that would be another factor that might help to characterize real events vs false alarm. In that case, better modelling might be highly desirable, but as a rough idea:
Imagine your three lists as three parallel lines in three-dimensional space. They are arranged with equal distance to each other, i.e. there's a fourth parallel line in their "middle" that has the same distance to each one of them.
Now, imagine your timestamps are points on the respective lines.
For every point on your magical middle line, you draw a sphere around that point and count the time stamps in that sphere. Setting an appropriate sphere radius and threshold allows for good detection.
If you have a mathematical model for how observed timestamps relate to real phenomena, that can change the shape of your sphere, and make it "fuzzy", i.e. weighing points further away less than nearby.
In the end, you'd have a function that says "for each time on my middle line, this is the likelihood of a real event having happened"; you can then build a maximum likelihood detector.

¹ hint: This might be a bad stereotype, but you sound like a matlab user. If there's one easy thing matlab is slow at, it's doing actual m-code loops and removing heads from lists; it's always preferable to just work at the end of a contiguous vector of data, since all you need to do to "remove" an element is change your own notion of length.


Finding the amplidude of a single frequency in signal I'm not experienced whatsoever in signal processing. I want to isolate a single frequency (83 Hz) out of a noisy signal. Goertzel works just fine but was wondering if there was any other, maybe even (computationally) simpler technique to do that. I only want the amplitude of that single frequency.

Goertzel in fact is the matched filter for a single frequency – so, it's, in the presence of uncorrelated noise, the estimator that gives the best estimation variance under fixed observation.

But: for Goertzel to work, you need to use a number of samples that is

an integer and
a multiple of $f_\text{sample}/83\,\text{Hz}$.

Unless your sampling rate hence is a small multiple of 83 Hz, that's usually a bad deal.

You can remedy that with a number of things: often, a narrowband (often: even resonant) filter can filter out a frequency range, but you'll need to deal with the skirt of that filter. That skirting will be different than that of Goertzel's algorithm (unless you happen to have also windowed the signal before applying the Goertzel, and now use the same Window as filter design prototype; sounds unlikely).

If you know more about the signal than you're telling us (it's result of AR processes, so you can apply Yule-Walker to get an appropriate model of what's happening, for example. Or, you know the signal is the sum of white noise and a known number of oscillations – then, autocorrelation subspace methods such as MUSIC are very efficient), then you might have methods that are very exact, but use even fewer samples than minimal Goertzel.

Anyway, the beauty of Goertzel really lies in the fact that it's pretty efficient to begin with. Honestly, if your frequency of interest is 83 Hz, and you still can't detect it in realtime even on the weakest of microcontrollers, then you're doing something wrong, with 99.5% probability¹. You really might want to ask a different question on here, describing what you overall want to do (where does your signal come from? What's at 83 Hz? Why do you want to determine the amplitue of that? How exactly do you know that frequency? What's the amplitude estimate reliability that you really need? So many things to talk about...).

¹ This number, like 22.36% of all study results, is totally made up.


Confusion about the impact of flat fading and frequency selective fading on OFDM OFDM divides the signal in frequency domain into orthogonal and relatively narrow subcarriers where fading can be considered to be flat
Flat fading is not associated with delay spread and is corrected in frequency domain

Delay spread caused by multipath is associated with frequency selective fading and in OFDM it is absorbed by the cyclic prefix which isolates consecutive OFDM symbols

How would delay spread show up in the frequency response of the channel? would it show up like "nonflat" fading or a notch within the relatively small bandwidth of some subcarriers?

i.e. it would be wrong in this case to consider the affected subcarriers still experiencing flat fading?

Because of wide-band nature of OFDM signals they are so prone and sensitive to multi-path delayed channel. The problem is that at the receiver side you will receive several copies of original transmitted signal, each with a different magnitude and delay. This "delays" are so destructive at the receiver side. The mathematics behind this: imagine the channel with multi-path problem. so u can model this channel as a linear-system with an "impulse response". this impulse response consists of several impulses with different magnitude and delays. so, your sent OFDM signal passes through this channel. the result can be guessed. at the output this system, the OFDM signal is convolved with the impulse response of the channel (system) in time domain. In frequency, it is like as the freq-domain OFDM signal is multiplied by the frequency-response of the channel. That's why we face a distorted received signal. At the receiver, when you apply FFT on the received signal you have to compensate that errors. This is called channel estimation (means estimating the freq-response of the channel). Then you should correct the signal. This is so crucial because in OFDM we usually have QAM modulation & demodulation. and in QAM the magnitude of signal is important which is already distorted by the channel.


Very basic question about filters in signal reconstruction I am studying digital signal processing and I would like to ask a question about sampled signal -> original signal reconstruction filters.

Most textbooks use a box filter or a tent filter to show how a discrete-time sampled signal can be reconstructed to (an approximation in this case) of the original signal. However my question is:

Given a discrete signal result from sampling a signal with a sampling interval "Ts", the spectrum will consist of infinite replicas spaced at 1/Ts and scaled by 1/Ts.  Now I assume 1/Ts is more than double the original signal maximum frequency, so there is no aliasing.

Our intent with the box (rect) or tent filter is try to "cut off" the spectrum replica in the middle with convolution. What most textbooks do not explain, is what size we should choose the rect or tent to be. Most examples I've seen use a box filter which has a base of Ts, or a tent filter with a base of 2*Ts. Why the base of the tent filter is double of that of the box filter? Is this an arbitrary choice? Looking at the fourier transform of a tent function (squared sinc) compared to that of a rect function (sinc) we shouldn't need to double the base to get the cutoff frequency near the (ideal) 1/2*Ts. With a tent with base 2*Ts  we actually hit 0 amplitude at 1/4*Ts. Aren't we cutting out more frequency than we should?

Neither a filter with a rectangular impulse response, nor one with a triangular impulse response is a very good interpolation / reconstruction filter. If you use a filter with a rectangular impulse response, the resulting output is just a piecewise constant signal. Since you want the constant pieces to have exactly the duration between two samples, the width of the rectangle must equal $T_s$. You could call that a zeroth order interpolation.

A first-order or linear interpolation between the sample values is achieved with a triangular filter kernel. If you look at it in the time domain you can see that the triangle with its apex aligned with a sampling instant $k$ must extend from the previous sample $k-1$ to the next sample $k+1$ in order to provide linear interpolation when all those triangles are summed up. So for linear interpolation, the base of the triangle must have a width of $2T_s$.

Of course, that triangular function is the convolution of two rectangles of width $T_s$, so its Fourier transform is the square of the Fourier transform of the rectangular pulse.


How do I interpret these autocorrelograms? Hello signal processing wizards.  I am young and feeble, and would like some help interpreting my spike-train autocorrelograms of six isolated neurons.  These histograms were created using the code pasted below (hist(k)).  Rather than creating an autocorrelation of the actual spike trains, this code correlates interspike intervals (see the following post for another example:

https://stats.stackexchange.com/questions/148865/how-to-interpret-this-autocorrelogram-graph

I am glad to have gotten the code working, but now I am struggling to interpret what the histogram actually means.  I am aware that the Y-axis shows how strongly the signal is correlated with itself at a particular lag (the X-axis).  However, I am not certain how I can relate these findings to the firing pattern of the neuron.  For instance, if I see a peak at say, ~ 750 ms, does this mean that with a lag time of 750 ms, the firing pattern is more constant/predictable, or does it simply mean that I have more inter-spike intervals of 750 ms?  If someone could clear up my confusion, that would be much appreciated.  On a related note, if anyone has any tips on how to choose an appropriate bin size and window, that would be awesome.

P.S. The X axis ranges from -1000 ms to 1000 ms.  Bin size is currently set to 50 ms.

time_of_spike= sparse([],[],[],1,tmax/dt);

clus_ms_Mult_10 = round((clusterTimes*10));%cluster time-stamps were in ms

while index < length(clusterTimes);
    time_of_spike(clus_ms_Mult_10(index))=1;
    index=index+1;
end

% % %finding the interspike intervals matrix
I=find(time_of_spike)/10;
i=size(I,2);
for m=1:i
    for n=1:i
        ISI(m,n)=(I(m)-I(n));
    end
end

%calculating the autocorrelation
b=1;
clear k;
for i= MIN:1:MAX
    k(b)=size(find(ISI(:,:)>=i & ISI(:,:)<=i+50),1);%if i want to change "bin
                                                    %size", I can change i+1 to i+largerInt
    b=b+1;
end

To answer your specific question:

  For instance, if I see a peak at say, ~ 750 ms, does this mean that with a lag time of 750 ms, the firing pattern is more constant/predictable, or does it simply mean that I have more inter-spike intervals of 750 ms?

You are calculating a histogram of inter-spike distances, binned at 50 (in the code as written). Obviously the spikes are only zero time instants away from themselves which is why the high peak at zero.

If there is a peak in the histogram at 750ms, it does not say anything about repeated inter-spike intervals of 750ms. That is, just because there is a second spike 750ms after the first spike, it doesn't necessarily mean there will be a third spike 750ms after the second.

It just says that, over time, there is more likely to be a spike 750ms after any given spike.


matlab's FFT is wrong? don't think so i'm trying to get an FFT of a signal

my FFT has 1024 samples

i added two sin functions one of freq 86,13Hz (2*44100/1024) and one of freq 129,19Hz (3*44100/1024)

the WAV looks like that (for 1 period) :

i then use this code to inspect the FFT

FFTSize = 1024;

Y = wavread('2sins.wav');

current_window = zeros(FFTSize, 1);

for k = 1:length(Y)/(FFTSize)

   for m = 1:FFTSize
       current_window(m, 1) = Y((k-1)*(FFTSize)+m, 1);
   end;

   plot(current_window)

   a = abs(Ydft);
   p = angle(Ydft);

  end;

the bins 2 and 3 have the maximum amplitude of 256 which is right as i used two sin of frequency 2 and 3

but when i look at the phases they are -pi/2 which is weird as they should be 0 no ?

here are the phases :

so you see they are -pi/2, why is that ?

and here is the WAV file i used so you can try yourself:
http://www.khaelis.com/tmp/2sins.zip

thanks for help

Jeff

Matlab is right, as usually. The DFT of a sine is imaginary. If you do the same with the sum of two cosines, you would get a real-valued result.


How to calculate the impulse response of an RC circuit using time-domain method I have an $RC$ circuit and I calculate the differential equation as
$$RC \frac{dV_c}{ dt} + V_c = V_g$$  that is the same of
$$RCy'(t) + y(t) = x'(t).$$
 If I assume $x(t) = u(t)$ I have that the equation is  $$RCg’(t) + g(t) = 0$$  and the solution is  $$g(t) = K_1 e^{-t/RC} + K_0$$  now for have the impulse response I have to calculate this for $t > 0$, $t < 0$ and $ t= 0$ but I have no idea how to do. Can you please help me ?

Dynamic electric circuits involving linear time-invariant resistors, capacitors, and inductors are described by linear constant coefficient differential equations (LCCDE). Mathematical solution of such LCCDE requires some physical (electrical circuit theoretic) insight too. This is especially true for solving circuits under impulse functions (such as finding impulse responses). Note, that the dirac impulse $\delta(t)$ is a non-differentiable, non-continuous function; therefore it cannot be the subject of a continuous-differential equation. Nevertheless, a procedure as I present below, is developed to solve LCCDEs involving such improper impulse functions.
Given the circuit assuming zero initial conditions, derive its DE by applying a KCL (Kirchhoff's current law) at the output ($y(t)$) node :
$$ \frac{y(t)}{R} = C \frac{d}{dt} \left( x(t) - y(t) \right) $$
$$ y(t) = RC x'(t) - RCy'(t) $$
$$ \boxed{ y'(t) +\frac{1}{RC} y = x'(t) } \tag{1} $$
Impulse reponse $h(t)$ is the solution of Eq.(1) when the input is an impulse: $$ x(t) = \delta(t) \implies y(t) = h(t)  $$
Let's apply a two-stage procedure to find $h(t)$. First, modify RHS of Eq.(1) by replacing $x'(t)$ with $x(t)$ :
$$ y'(t) +\frac{1}{RC} y = x(t)  \tag{2}$$
and denote solution of Eq.(2) as $h_0(t)$ when $x(t) = \delta(t)$. Then, by the second stage, the overall impulse response $h(t)$ of Eq.(1) will be:
$$\boxed{ h(t) = h_0'(t) .} \tag{3}$$
Solution of first stage, Eq.(2):
From the characteristic equation : $s + \frac{1}{RC} = 0 \implies s = - \frac{1}{RC}$.
The homogeneous solution is :
$$y_h(t) = K e^{-t/RC} u(t) \tag{4}$$
The particular solution $y_p(t) = 0$, as the input is an impulse; i.e., $$\lim_{t \to \infty} x(t) = 0 \implies y_p(t) = 0$$
And the complete solution:
$$h_0(t) = y_h(t) + y_p(t) =  K e^{-t/RC} u(t) $$
To find the unknown $K$, we can use $y(0^+) = h_0(0^+)$ ; and from the circuit physics we find $y(0^+) = h_0(0^+) = 1$.
Then $K = 1$, and the total solution of Eq.(2) is
$$ \boxed{ h_0(t) = e^{-t/RC} u(t) } \tag{5}$$
$u(t)$ function is used to indicate that the output is causal; i.e., zero for $t < 0$. (a realizable physical circuit will be a causal system.)
Finally, the impulse response $h(t)$ is found from Eq.(3), by the second stage, as:
$$ h(t) =  h_0(t)' = \left( e^{-t/RC} u(t) \right)' $$
$$ \boxed{ h(t) = \delta(t) - \frac{1}{RC} e^{-t/RC} u(t) } \tag{6}$$
An intuitive physical explanation of the circuit behaviour under impulse excitation is as follows. Due to the physical capacitor's continuity of voltage  property, the initial impulse (physical spike) voltage applied at the input, appears right at the output as a result of KVL applied around the loop at $t=0$. This impulse voltage on the output resistor $R$ produces an impulse current $i(t)= \delta(t)/R$ at time $t=0$ looping clock-wise around the circuit. This impulse current passes through the capacitor and jumps its voltage from $V(0^-) = 0$ to $V(0^+) = 1/RC$. Then the input node is shorted as the impulse is gone at $t=0^+$. And the simple first-order capacitive discharge circuit is obtained. Note that the discharge current loops counter clock wise due to capacitor voltage polarity. Hence the output voltage becomes negative for $t>0^+$.


The model of noise in Image denoising problem I am studying the image denoising problem. The paper I studied uses Gaussian Noise to model the noise appearing on the image. E.e., If we say $u_g$ is a clean image, then we model the noisy image $u_0$ so that $u_0 = u_g+n$ where $n$ has Gaussian distribution with mean $0$ and variance $c$.

Q1: I am wondering why are we using Gaussian distribution to model the noise? Any theory behind it?

Q2: Can we use stochastic processes to model the noise? If so, is there any existing literature?

Theory behind it: the additive Gaussian noise  can be related to (Johnson–Nyquist) thermal noise in standard images. For the  denoising of more generic multivariate images, chapter 2.2.2 provides some references on why global noise present in acquired data can be realistically modelled by an additive zero-mean spatially white Gaussian noise.

The shot noise can be modeled as a Poisson process, that behaves, for large numbers, like a Gaussian noise.

Other types of noises are observable in images: salt-and-pepper, quantization (more uniform), grain. In some image modalities, noise can be not stationary, multiplicative. They can be modeled by stochastic  processes in general.

Gaussian noise is also widely used because it enjoys relatively tractable theories, related to least-square estimation. A very practical reason.

Note, however, that most images have 8-bit channells, with integer values in $[0,\ldots,255]$. hence, adding a Gaussian noise yields non-integer negative or values above 255, which is not  fully realistic: natural sampled  noise is  likely to be quantized in the  admissible values.


Filtering White Noise to Generate Harmonic Signal I have studied convolutions and filters a long time ago. Today, I am trying to review the basics using some notes of mine, but I am finding difficult to solve easy problems. Since I don't have solutions for these problems, I believe that it could be helpful for me if you could show me how to solve one of them.

Let's say that I want to define a filter $b_\alpha(L)$ (where $L$ is the lag operator) that applied to a white noise $x(t)$ approximates a $\cos(\alpha t)$ function (where $\alpha$ is given). How should I proceed to estimate this filter?

If you filter white noise with a narrow band filter that passes $e^{j\alpha t}$ and $e^{-j\alpha t}$ you will get a sinusoidal output at mean frequency $\alpha$ with an amplitude and phase that varies at the maximum rate according to the bandwidth of your filter.  The narrower the bandwidth, the lower the power (since the power for white noise is a distribution evenly spread over all frequencies), but the closer you will approximate a pure cosine signal at a random initial offset in phase and random initial amplitude (which will then still proceed to vary, although relatively  slowly compared to a winder band filter, in amplitude and phase).  If you could do an infinitely narrow filter (which you can't), you would also get no signal.


Techniques to deriving DTFTs Are there general techniques to derive DTFTs? Given a bandlimited function $x(t)$, how do I find

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

Generally, it is easier to derive the continuous transform (never mind the constants):

$$X(f)=\int_{-\infty}^{\infty}x(t)e^{-i \omega t}\, dt$$

because we have a wealth of integration theory to fall back on. Ok, one can also calculate $X(\omega)$ from $X(f)$ by adding shifted copies, but we're again back to summing an infinite sequence. Any ideas on how to approach DTFT derivations in a smarter way? Or is sweating out the tiring summation the only way? I don't have any particular function in mind, but if one has to be chosen as an example, I'd pick $x(t) = \text{sech}(a t)$. The continuous FT of this is known and can be found in the table on Wikipedia (#208).

As an example, the rules of series can work. A typical example is the $x[n] = a^n u[n]$ where $|a| < 1$. For instance, say we want to find the DTFT of the signal $x[n] = 0.9^n u[n]$. Then,

$$X(e^{j\omega}) = \sum_{n = -\infty}^{\infty} 0.9^n u[n] e^{-j\omega n} \\
= \sum_{n = 0}^{\infty} (0.9 e^{-j\omega})^n$$

This is where the series work (geometric power series):

$$\sum_{n = 0}^{\infty} a^n = \frac{1}{1-a}$$

Then our DTFT is:

$$X(e^{j\omega}) = \frac{1}{1 - 0.9 e^{-j\omega}}$$


Video Compression - When would an average of the previous and next I or P frame be useful I'm doing practice exam questions for my multimedia systems exam one tuesday and a topic which seems to come up frequently is the appropriate use of encoding B-frames based on:

1)  the previous I or P frame

2)  the next I or P frame

3)  an average of the previous I or P frame

I'm having a bit of trouble thinking up decent examples of these, and I'm not entirely sure when an average would be most appropriate.  Would someone mind checking my examples or maybe helping me think of some more concise or realistic ones?

1)  A dropped ball about to hit the floor.  The previous frame the ball is the same shape but in the next frame it gets deformed because of the impact of hitting the ground, so the previous frame is better.

2)  A man running along the pavement and a car moves between him and the camera temporarily in the previous frame.  Due to object occlusion the image of the man in the previous frame is not complete so is not a good predictive frame, but the next frame the image of the man is complete so it is better to predict from this.

3)  I'm not sure at all what to put on this one.  I was thinking about a car moving with constant velocity across the screen because then it's position in the target would be the average of the previous and next -  but this seems quite contrived!

Anyway, any help would be appreciated so if you have ANY thoughts let me know!

Before elaborating anything, i would urge you to go through my answer here that simplifies some of the common confusions prevailing about motion vectors and estimations.

Your essential question is - when or how does encoder decides to put B frame or bidirectional motion vectors. I will break your questions into few as follows:

1. When does encoder decides to use B frames?

Almost always the number of B frames between I or P is always fixed. This is usually 2 or 3 pictures. The distance between I frame (called GoP size) may keep varying, but number of B frames will always remain fixed. From the point of view of MPEG, this is not really a restriction, but as one of your other question discusses, if you keep changing B frames, there will be serious problems of rhythms between the transmission and receptions.

Also, even before estimating whether a given frame will be a good choice to act as a B frame or a new I or P frame, you have to actually decide the forward reference and past reference for from where the prediction. Iterating over whether given frame will be perfect to work as P against frame already gone or should you choose some distance frame as P and current as B is too much overhead to get encoding done, and in the process you are fluctuating delay heavily.

So in essence, typically encoders decides number of B frames as 1,2 or maximum 3 and corresponding pictures are always fixed to be B frames.

2. So then what if i have a huge black/uncorrelated picture between last I and future P?

Ok - now that you have chosen to see that given frame is B, but almost nothing matches of that frame between either past or future; or sometimes only either side of the frame? This doesn't really put encoder in trouble though.

In MPEG, every frame is divided in terms of macro blocks. Each macroblcok has motion prediction from forward or reference blocks from respective past or future pictures.

I frame can ONLY have Intra blocks - which has no motion compensation, so it doesn't need any other frame ready for decoding itself.
P frame can have INTRA as well as FORWARD-PREDICTED macro-blocks. So even in case of P frame, when current block doesn't match with anyone in the past frame, P frame can simply use INTRA block.
B frame allows all type of macroblock including INTRA (as in I), FORWARD-PREDICTED mb (where motion vector depends only on past picture) , BACKWARD-PREDICTED mb (where motion vector depends only on future picture) and BIDIRECTIONAL macroblock which has two motion vectors from past and future and they are averaged out to form current block.

Hence, in case, the B picture or part of it happens to be very independent it is quite possible to put INTRA blocks for those MBs, and if it matches only either side of the picture encoder can use only FORWARD or BACKWARD prediction either way.

3. So does all MB's of a given frame have to be of same type?

Not at all, suppose a new small object just entered in a frame which is encoded as B frames, only part of the frame needs INTRA Macroblocks but majority of others can predict their motion from other blocks.

Understand, that all MPEG cares about is, if encoder can find a block which is almost as good as a current one, use to minimize the error, if not bad luck and encode this as a fresh new information. Even this ability to find itself varies from encoder to encoder.

4. So how does encoder make so many choices?

Usually, for the B pictures, encoder starts a usual business by searching the most appropriate blocks from past and followed by future picture.

a. if both frame produce some reasonable motion vector, this is most preferred to minimize the energy and book the block as BIDIRECTIONAL MB.

b. if either past frame or future frame doesn't produce motion vector - i.e. the minimum energy (SAD) is higher than a useful number, then correspondingly only 1 motion vector is encoded making MB FORWARD or BACKWARD prediction only.

c. if none of the reference frames have good enough motion vector - the MB is encoded as INTRA.

5. So why don't I simply encode everything as INTRA block and save from all the hassel?

You can absolutely do that. It will still be perfectly legal MPEG stream. Just that INTRA coding will take more symbols to encode and hence it will have to be quantized badly to fit under a given bit budget. Putting everything as INTRA has significant loss of efficiency because typically video always has huge temporal redundancy that INTRA blocks don't exploit; but they do give random access freedom.


Find signal power in a band by integrating the PSD in a frequency band (and zeroing the out-of-band DFT bins) I know that directly zero DFT bins outside a frequency band has the side effect of introducing ringing, as this post says Why is it a bad idea to filter by zeroing out FFT bins?. But what about calculating the power by summing the PSD in a frequency band. In my opinion, this is similar to first setting the out-band PSD amplitudes to zeros (which may cause ringing in time domain) and then summing the remainning nonzero amplitudes.
I guess it is different from first applying a bandpass filter to the time domain signal(which has no ringing) and calculating the power by averaging the squared time domain samples. Can you point out how much error the frequency domain method has? or maybe these two methods are equivalent?

For a DFT that was computed with only a rectangular window (no further windowing beyond the selection), the power in the "out of band" bins contains signal energy from "in-band", so by zeroing those out, you are not including the power that should be in the computation. The challenge however is that if there is energy out of band, it will similarly have components within your band of interest! Here we see that the real flaw is trying to comput the power within a bandwidth without windowing the signal in the time domain first to reduce the effects of sidelobes in the DFT.
More details on doing this accurately, and specifically accounting for the power loss due to the windowing, and the change in the equivalent noise bandwidth of each bin, please see this post here:
How can I get the power of a specific frequency band after FFT?
and
How to calculate resolution of DFT with Hamming/Hann window?


What is the meaning of the sum of kernel components in digital image filtering? When analyzing an image signal as a series of Dirac impulses in the continuous spatial domain, I applied continuous low pass filters (to reconstruct or smooth it) with area equal to 1 to preserve the low frequency amplitudes in the signal spectrum. When using a kernel (which, correct me if I'm wrong, is just the "sampled" discrete version of the original continuous filter in the spatial domain) I was wondering what is the meaning of the sum of the kernel components, and when I should care?

I often see examples with kernels with components which sum to one, stating they preserves the original image brightness. But if the original continuous signal integrates to one, it does not mean the discrete version will sum to one, so this left me a bit confused.

Considering an LTI filter impulse response (aka filter kernel) $h[n_1,n_2]$, then its frequency response will be

$$ H(\omega_1,\omega_2) = \sum_{n_1} \sum_{n_2} h[n_1,n_2] e^{-j \omega_1 n_1}e^{-j \omega_2 n_2} $$

When an image is filtered, its average brightness is mostly affected by the DC response of $h[n_1,n_2]$, which is given by $H(0,0)$. When you want the brightness not to change after filtering, then you set $H(0,0) = 1$, which means that

$$ H(0,0) = \sum_{n_1} \sum_{n_2} h[n_1,n_2] e^{-j 0 n_1}e^{-j 0n_2} = \sum_{n_1} \sum_{n_2} h[n_1,n_2] = 1 $$

Which shows that the sum of impulse response (aka filter kernel) samples is one.


What Are The Semantics Of Wav-File Sample Values? The title is the question.  To make this concrete, assume I have a standard 16-bit single channel wav file, with (although it doesn't matter) a 16000 Hz sampling rate.

Each sample magnitude will be a 16-bit value.

Is there a specified semantics or interpretation of these values, not relative to any physical value, but relative to the numerical systems of modern computers?  I.e., are those numbers to be interpreted as 16-bit twos-complement, or scaled into single precision floats between -1 and 1, etc?

At the highest level, I recognize this is a little pointless:  They're arbitrary scale.

But at a practical level, in the python world, librosa does one thing, scipy does another, and converting those different things into a spectrogram or analyzing with a vocoder can yield different results.

For skeptics, invite you to try the following:

import scipy
import librosa.display

np_eps        = np.finfo('f').eps
np_max        = np.finfo('f').max

x, sr = librosa.load(file, sr=16000)

plt.subplot(2, 2, 1)
D = librosa.amplitude_to_db(np.abs(librosa.stft(x)), ref=1.0)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')

plt.subplot(2, 2, 2)
plt.plot(x)

sr, y = scipy.io.wavfile.read(file)
y = y.astype(float32)
plt.subplot(2, 2, 3)
D = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=1.0)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')

plt.subplot(2, 2, 4)
plt.plot(y)

Do not neglect the extra argument to librosa, forcing the sampling rate; in my case, the file I am using is from the Arctic corpus and is sampled at 16 KHz.  Do not neglect the colorbars, which will show you the very obvious shift of 90 dB from one to the other.  This should not be surprising, since 20 * log10 (32768) is approximately 90.3.

I assure you, gentle readers, this will have consequences further down the pipeline of speech vocoders, or neural networks (where normalization techniques might want you to reference to numpy epsilon to avoid negative numbers) or both.

I don't know much about this semantics? of WAV files but their numerical format is the following. (assuming mono)

Given a recording with 8-bit per sample precision, then those samples are unsigned integers taking values between $0$ and $255$. Due to being unsigned, to represent negative values,  there is a bias of $128$, and the sample values are actually interpreted to be between $-128$ and $127$. In this case a value of $128$ represents silence, maximum value $255$ corresponds to an input full scale positive audio waveform and $0$ represents neagitve full scale input waveform.
For $16$-bit WAV files, you have $16$-bit signed integers that range from $-32768$ to $32767$ and $0$ represents silence. Again, $32767$ represents full scale audio input without scaling.

Note that here I refer to a microphone recording through an ADC system which converts analog audio input in the standard range of about $1$ V rms into WAV samples. Please check this 1-V rms for being something higher, about 1.5 Volts or else depending on particular implementation or some specific audio standard. In any case, the analog circuitry is responsible for transmitting a full scale audio to the ADC whan a 1V (or 1.5V) rms wave is present at the input and then the conversion is performed according to the above setting when WAV recording is performed.


Anyone explain to me this video? I was watching a
video
in time 24:48
I would like to know where you got the value

.9 (1.14z + .941) and 1.0232 + .757

Does anyone explain how he got those numbers?


Calculating the energy of trapezoid wave Here is the trapezoid wave,although i know we can use $\int^{4}_0 |x(t)|^2dt$ to calculate its power,but i want to divide this trapezoid to right triangle$(t=0$~$2)$,and a square ($(t=2$~$4)$).

I know the the triangle effective voltage is $V_{eff}=\frac{V}{\sqrt{3}}$,and i want to use $P=\frac{V_{eff}^2}{R}$ to calculate the power ,  but i am not sure about which  triangle is correct,i mean,for which triangle(a or b or c),i can use $\frac{V}{\sqrt{3}}$ to calculate its power,or are these three wave's power the same,but different in energy?

it seems that there power are the same,but there energy are not the same ,because energy=power $\times$ time,and time are different,so their energy are not the same too


"chirp" with arbitrary period Say you have a linear chirp, which is a bit like a sinusoid with a gradually increasing period, but instead of the linearly increasing period, could you pick an arbitrary value, like the red line in the image below?

How would you write the definition of this? Say that the period variable is called AP (arbitrary period).

For the linear chirp I refer to wikipedia.

If I understand correctly, you are trying to create a chirp an arbitrary instantaneous frequency, $f(t)$.
In that case, the phase of the chirp is given by $\varphi(t) = 2 \pi \int f(t)dt$ and your signal in time domain is given by $s(t) = \sin(\varphi(t))$.
For example, in case of a linear chirp  you mentioned, $f(t) = f_0 + \beta t$ (where $\beta$ is the sweep rate) and therefore $s(t) = \sin \left(2 \pi (f_0 t+\frac{1}{2}\beta t^2) \right)$.


Ultrasound pulses and reflection problems I am building a whiteboard to track the position of its pen. I have n ultrasound receivers placed on the periphery of the whiteboard, and an ultrasound emitter in the pen. The pen emits pulses the receivers detect. A microprocessor gathers the arrival times of the pulses for each receiver, and with the time differences of arrival (TDOA) it makes an estimate of the position of the pen following a multilateration algorithm.

Ultrasounds reflect on surfaces, which can cause interference. This is a restriction on the frequency of the pulses. At the moment, I pulse the pen at 10Hz, which is enough for the reflections to die off and not interfere. Unfortunately, 10Hz is not fast enough a sampling for my purposes. Ideally the pen would pulse at 100Hz. (The pen emitter is connected to a microcontroller, so I can control the shape and frequency of the pulse.)

What tricks can I use to tackle the reflection issue? What are some standard filtering techniques? Could using different pulse patterns in a cycle help filter out the reflections?

It would appear that your problem lends itself very nicely to using a CDMA scheme.

Let starts with some properties of (DSSS) CDMA. (Direct Sequence Spread Spectrum, Code-Division Multiple Access). Its a mouthful, but it is really easy to implement.

In CDMA, your pulse (at baseband) is actually made up of many concatenated 'chips' as they are called. The chips are just 1s or -1s, of a fixed duration. For example, your chipping sequence might be [1 -1 1 -1 -1 -1 1]. You would use this chipping sequence to modulate your carrier.

However, you cannot just make up your chipping code. What you want to do is use chipping codes that have the very nice property, that their autocorrelation function is a delta function like so:

(Equivalently, their power spectral density is white). For example, you can look into using Barker Sequences as your chipping code, (usually used in radar), or you can also look at using Gold Codes. Practically speaking however, this means that you get the maximum correlation score in your receiver, ONLY when the receivers' code, exactly lines up with the transmitted code, and zero otherwise.

How does this help you? In your receiver, you would be running a correlator continuously. The correlator would be performing a running dot-product of its own local code, with whatever is received. Now imagine that you receive a transmitted waveform from your pen, and a second waveform from a reflection. As your receivers' correlator runs, it will give a peak when its own codeword exactly alligns with your code from the pen. This will cause your detector to 'lock' onto that specific delay value. Now, here is where you reap the benefits of a near-delta autocorrelation function of your code: The reflected signal will also be present, and will also have its dot product taken with the receivers' locked code, but it will give zero, or near zero score, since it is orthoginal or near-orthogonal to the delayed code that your receiver has already locked onto.

In contrast, if you had send out a un-coded carrier pulse, you would be at the mercy of constructive or destructive interference throwing off when exactly your pulse peaked at the detector level of your receiver, and thus get erroneous TDOAs.


Machine Vision Techniques for Sign Language Recognition and Interpretation(using Matlab) I'm currently working on applying Machine Vision methods to sign language interpretation and I'm really overwhelmed with the body of knowledge out there. I've been reading up on object recognition from various resources (Peter Corke's Robotics, Vision and Control, Richard Szeliski's Computer Vision, Algorithms and Application), but finding it hard to decide on what methods to go with, hence my question:
For simplicity's sake, let's assume I just want to process an image of a hand holding up 1 finger, get my program to process the image and out "one". Do same for 2...up till five. (I know this is really simple, but will give me the direction I need before going on to more complex things).
(1)What method do I need to store the original image of the hand holding up 1 finger to put in my database, which other images will be compared to?
(2)What method do I to compare an inptu to what's stored in my database (Please, give me the most simple methods you know of).

Thanks in advance.

I am afraid that the question is trying to probe a very wide range of issues to receive a truly meaningful answer.

In general, this would be a classification problem. As such, you would need to decide on a number of features which would reflect the key changing characteristics of your images and a classifier which, after being trained, would be able to discriminate between your different gestures. For more information please see: http://en.wikipedia.org/wiki/Statistical_classification

The issues of how many training images you need per gesture and which features and classifier, all depend on your environment and specifics of your application.

The simplest thing you can do is template matching, where you would simply compare an acquired image of a gesture with stored images of gestures. For more information please see: http://en.wikipedia.org/wiki/Template_matching

But of course, this approach is a tiny step forward. What about changing backgrounds, different pairs of hands, scaling of the objects because of different distances to the camera, slight rotations, occlusions... ? And then, once the static recognition is done, you will have to consider the version that works on uninterrupted video.

There are things you can do to deal with these issues but for robust recognition, it is likely that you will need to move to 3D data representing the movement of hands in space. (For example: http://www.sciencedirect.com/science/article/pii/S1077314210001748 ).

You can acquire 3D data in a number of ways, e.g. virtual reality gloves (please see: http://dx.doi.org/10.1109/IJCNN.1999.832699), devices such as kinect (please see: http://dx.doi.org/10.1109/ICSESS.2012.6269439) and plain simple accelerometers (please see: http://dx.doi.org/10.1109/JSEN.2011.2166953).

The simplest method to start with 3D data is the one based on accelerometers. In this case you would be working with time-series data (e.g. timestamp, acceleration vector).

Happy to amend the response if you provide more specific data about your application.


Whitening Matched Filters for Equalizer I was looking at the answer Whitened Matched Filter. It mentions that in case of ISI channel,  we can use $g^∗(−t)$ as matched filter. For ISI channel, $r(t)=c(t)∗s(t)+n(t)$ where $r(t)$ is the received signal, $c(t)$ is the channel response, $s(t)$ transmit signal and $n(t)$ is the noise signal. The receiver structure for ISI channel consists of matched filter $g^*(-t) = c^*(-t) * s^*(-t)$ followed by a sampler (A/D).

My question is, that if this receiver structure practically implemented or is it used for analysis only? Since in a real receiver, the matched filtering is implemented at the front-end does not include the channel response but only signal response ($g^*(-t) = s^*(-t)$). This is because we don't know the channel response before the sampler and we are not adapting analog filters to match to the channel response.
How does a practical receiver structure looks like with adaptive filters after the sampler on the digital side? What am I missing here?

Thanks
strong text

You are right that in scenarios where the channel is unknown and/or time-varying, we cannot implement a fixed matched filter. The first filter in the signal chain is usually an analog band pass filter eliminating out-of-band  noise. Then after A/D conversion one would use an adaptive equalizer that performs part of the function of the matched filter, but that will usually also be optimized for further noise suppression. This optimality criterion for equalizers that combines minimization of intersymbol-interference and noise power is commonly referred to as minimum mean-square error (MSE) criterion.


Locally Adaptive Thresholding - How to Detect if Image Block is Uniform I need to implement a locally adaptive threshold algorithm on a grayscale image in matlab. I use a sliding window over the image and calculate the threshold for the block using Otsu's method, which works well for blocks with non-uniform intensities.

However, when I apply this method to blocks with uniform intensities, I get very busy (junk) blocks. I want to set these blocks to either black or white in the result. I can decide whether to use black and white by comparing the global threshold value to the mean of the block (is this a good approach?).

My main question is: how to decide whether a block is uniform or not, or rather, how do I know which blocks to threshold and which to completely set to black or white?

If you calculate the standard deviation of the blocks (from the average) it could serve as a good measure to the smoothness of a block (assuming there is no an extraordinary value and the specific use of the "system")


What Are the Best Algorithms for Document Image Thresholding (Example Inside)? I'm trying to implement various binarization algorithms to the image shown:

Here's the code:

clc;
clear;
x=imread('n2.jpg');     %load original image

%   Now we resize the images so that computational work
  becomes easier later onwards for us.

size(x);
x=imresize(x,[500 800]);
figure;
imshow(x);
title('original image');

z=rgb2hsv(x);       %extract the value part of hsv plane
v=z(:,:,3);
v=imadjust(v);

%now we find the mean and standard deviation required for niblack and
%sauvola algorithms

m = mean(v(:))
s=std(v(:))
k=-.4;
value=m+ k*s;
temp=v;

%   implementing niblack thresholding algorithm:

for p=1:1:500
    for q=1:1:800
        pixel=temp(p,q);
        if(pixel>value)
            temp(p,q)=1;
        else
            temp(p,q)=0;
        end
    end
end
figure;
imshow(temp);
title('result by niblack');
k=kittlerMet(g);
figure;
imshow(k);
title('result by kittlerMet');

%   implementing sauvola thresholding algorithm:

val2=m*(1+.1*((s/128)-1));
t2=v;
for p=1:1:500
for q=1:1:800
    pixel=t2(p,q);
    if(pixel>value)
        t2(p,q)=1;
    else
        t2(p,q)=0;
    end
end

end

figure;
imshow(t2);
title('result by sauvola');

The results I obtained are as shown:

As you can see the resultant images are degraded at the darker spots.Could someone please suggest how to optimize my result??

Your image doesn't have uniform brightness,so you shouldn't work with a uniform threshold. You need an adaptive threshold. This can be implemented by preprocessing the image to make the brightness more uniform across the image (code written in Mathematica, you'll have to implement the Matlab version for yourself):

A simple way to make the brightness uniform is to remove the actual text from the image using a closing filter:

white = Closing[src, DiskMatrix[5]]

The filter size should be chosen larger than the font stroke width and smaller than the size of the stains you're trying to remove.

EDIT:
I was asked in the comments to explain what a closing operation does. It's a morphological dilation followed by a morphological erosion. The dilation essentially moves the structuring element at every position in the image, and picks the brightest pixel under the mask, thus :

removing dark structures smaller than the structuring element
shrinking larger dark structures by the size of the structuring element
enlarging bright structures

The erosion operation does the opposite (it picks the darkest pixel under inside the structuring element), so if you apply it on the dilated image:

the dark structures that were removed because they're smaller than the structuring element are still gone
the darker structures that were shrunk are enlarged again to their original size (though their shape will be smoother)
the bright structures are reduced to their original size

So the closing operation removes small dark objects with only minor changes to larger dark objects and bright objects.

Here's an example with different structuring element sizes:

As the size of the structuring element increases, more and more of the characters is removed. At radius=5, all of the characters are removed. If the radius is increased further, the smaller stains are removed, too:

Now you just divide the original image by this "white image" to get an image of (nearly) uniform brightness:

whiteAdjusted = Image[ImageData[src]/ImageData[white]*0.85]

This image can now be binarized with a constant threshold:

Binarize[whiteAdjusted , 0.6]


best robust feature detectors for detecting similarities among archeological signs I would like to implement edge detectors for creating a feature vector for archeological signs.
Which are best features, algorithms to retrieve? That could be loaded in a Support vector machine for classification.

If there are in the image only signs with a uniform background it's better converting an image in a binary image and calculate convexHull, solidity, eccentricity, elongation of the archeological signs.
Then, what i do? I insert these data in a vector and put every vector of every image in a Support Vector Machine for training?

How can i measure "distance", similarities among images?

There is a difference between edge detection and feature detection, which is usually a corner detection (e.g. Harris corner detector), blob detection or even region detection (e.g. MSER).

Once the features are detected in your image, they are described in some way. For point features (corners), there are many different descriptors.

I think the simplest yet powerful one is simply taking a patch around the detected feature location. The patch can be quite large and should be rotated in the direction of dominant orientation (for rotational invariance).

The feature vector can consist of just horizontal and vertical image gradients. For a 64 by 64 pixel patch, this would lead to 4096 dimensional vector. This is quite large, but can be effectively reduced by PCA (Principal Component Analysis). This is quite slow but provides good quality and robustness of the descriptor comparable to SIFT or SURF descriptors.

The SURF descriptor is based on Haar wavelet and is also quite easy to implement - it can produce smaller descriptor and hence PCA is not necessary.

You can take a look on this paper comparing various feature descriptors.


Given a causal LTI system with a transfer function, determine if the system is an all-pass filter We are given a transfer function of a causal LTI system and must determine if it is an all-pass filter:

$$ H(z) = \frac{1 + 4z^{-2}}{4 - z^{-2}} $$

To the best of my recollection to determine if a system is an all-pass filter, we must first find the poles and zeros and see if they are conjugate reciprocal pairs. So working on that and finding the poles and zeros, we find that...

the poles are at $z = \tfrac{1}{2}, -\tfrac{1}{2}$

the zeros are at $z = 2j, -2j$

With that being said the poles and zeros are no conjugate reciprocal pairs. That being said, are there any exceptions or other criteria that would make a system an all-pass filter?

In an all-pass filter we should have $$H(z)H^*(1/z^*)=c^2$$
which means that the poles and zeros of the system $H^*(1/z^*)$ cancel the zeros and poles of $H(z)$.
As a result, an all-pass filter has conjugate-reciprocal pole-zero pairs. For this system:

$$H^*(1/z^*)=\frac{4z^2+1}{-z^2+4}$$
So
$$H(z)H^*(1/z^*)\neq c^2$$
You can also see that the poles and zeros are not conjugate and reciprocal: $\frac{1}{\left(\pm2j\right)^*}\neq\pm\frac{1}{2}$


About phase unwrapping I am not too familiar with phase-related problems so I need some suggestions or corrections in a problem that might be an old but easy one. Here is my problem: given a complex signal:
$$y(t) = e^{j\phi(t)}$$
where $\phi(t)$ is the instantaneous phase. Now we have the observations $y(t)$. Can we somehow reconstruct $\phi(t)$?

I understand the phase will be wrapped into $-\pi$ to $\pi$ when we calculate the phase from $y(t)$. For example, let $\phi(t) = 130t+150t^2 + 0.2, t\in[0,1]$ with a sampling frequency $f = 512$. When I try to reconstruct the phase from $y(t)$ using unwrap(angle(y(t))) in Matlab, i.e., the unwrapped phase, the result is shown in the following figure:

The second example is that $\phi(t) = 450t-100t^3 + 0.2, t\in[0,1]$. The result is

The reconstructed phases are different from the original ones. Does this mean that the phase cannot be uniquely reconstructed? Could anyone show me some explanations or suggestions? Thank you!

Edit: the code is as follows.

clear all
clc
N = 512;
t = (0:N-1)/N;
% phi = (130*t + 150*t.^2) + 0.2;
phi = (450*t - 100*t.^3) + 0.2;
s = exp(2*pi*1i*phi);
phase_rcst = unwrap(angle(s))/(2*pi);
figure
subplot(2,1,1)
plot(phi)
title('Phase of the original signal')
subplot(2,1,2)
plot(phase_rcst)
title('Unwrapped phase')

The problem comes from sampling the phase. For the given phase function, the changes in the phase from one value $t_i$ to the next one $t_{i+1}$ can become too large, so phase unwrapping will not work properly. If you simply removed the factor $2\pi$ in the exponent, the problem would be solved already (note that the factor shouldn't be there anyway according to definition of $y(t)$ in your first formula):

s = exp(1i*phi);

If you want to keep the factor $2\pi$, then use more sampling points.


Correct way to quantize accelerometer data with microsecond precision? I am analyzing large volumes of raw accelerometer data (~12 hours @ 200Hz) using Python. Some calculations (e.g. deriving Jerk) necessitate a time index with a regular interval. I have microsecond-precision time associated with my data but I can see variance in the +/- 80μs range. Which method of time quantization will help me achieve a regular time interval while introducing the least amount of error?

I'm a DSP novice so please excuse and correct me if I'm misusing any terminology

useless ramblings first

As an SDR person: 12h · 200 S/s = 12·3600s·200S/s = 8.64 MS is not a "large volume of data" (we usually deal with a multiple of that per second), so rest assured, most of your computation should probably be pretty fast to complete, which is good news!

Relating your Jitter to your sampling period

Now, with 200 Hz, your sample period is 5000 µs, so a 80 µs variance (something doesn't add up with the units here, it should be s², not s) would be an error of ca 1.6%, or -18dB. If you assume sufficient band-limitation, rule of thumb says that the resulting phase noise would be in the same order of magnitude, ie. -18dB.

It's generally a good idea to define first what noise level are you aiming for; otherwise you'll be optimizing something that is already good enough.

Of course, you could properly resample, with much less distortion, but frankly, accelerometer data doesn't sound like you necessarily have great SNR in the first place (but I might be totally wrong about that!).

Interpolator considerations

What you need, if you actually want to generate a perfectly periodically sampled signal, is an interpolator.

What confuses me is that especially jerk is a derivative measure, i.e. you'd be finding an interpolation of your signal more or less anyway; so, where's the problem in using a standard interpolation method (scipy comes with some) to get perfectly periodic samples, should you not choose to implement your own interpolating differentiator (hint: splines are of polynomial shape and that makes differentiation especially easy, so you can directly calculate a derivative from your data/time points, and be done with it, without first calculating a periodically-sampled signal representation).

Bandwidth considerations

If you know that your signal of interest is contained in frequencies significantly below 100 Hz, then consider this:

Assume your signal is a sum of individually delayed, scaled sines of different frequencies.

The point of maximum derivative, i.e. maximum steepness, is the point where our sampling time offset $\Delta t$ hurts the most.

A single sine has a maximum derivative of 1 (proof: the derivative of a sine is a cosine, which has a max of 1).

If each of your $N$ sines $s_i$ is scaled with an amplitude factor $a_i$ and has frequency $f_i$ (and phase==delay $\varphi_i$, which we do not care about right now), then its formula is

$$s_i = a_i \sin(2\pi f_i t + \varphi_i)$$

and its time derivative

$$\dot s_i = a_i\cdot 2\pi f_i \cdot \cos(2\pi f_i t + \varphi_i)$$

has a maximum of

$$\max \dot s_i = 2\pi a_i f_i = (\max s_i) \cdot 2\pi f_i \text.$$

In a linear approximation (and that's a prett good approximation for a range of less than 1.6%/2 = 0.8% of a cycle of a cosine around it's zero crossing; the non-scaled cosine is pretty much a 45° line there), the (defined-by-me) relative error (relation of error magnitude to amplitude of sine) you get is

$$\begin{align}
\max |e_i| &:= \max \left\lvert\frac{\dot s_i \cdot \Delta t}{a_i}\right\rvert\\
& = 2\pi f_i \Delta t\text,
\end{align}$$

which means that the lower the sine's frequency, the less you're sensitive to jitter.

Let's look at the case where you only care about frequency up to 50 Hz, or 5 Hz:

$$
\begin{align}
\max |e_i| &= 2\pi f_i \Delta t&f_i = 50\,\text{Hz}\\
&= 2\pi\, 5\cdot10^1\, \cdot\, 8\cdot10^{-5}\\
&= 2\pi\, 4\cdot10^{-3}\\
&\approx 2.5 \%\\[4em]
\max |e_l| &= 2\pi f_l \Delta t&f_l = 5\,\text{Hz}\\
&\approx 0.25 \%
\end{align}
$$

So, the lower your signal of interest's frequency is, the less you need to actually care about the jitter!


Problem with HSV and RGB color model I am a beginner in digital image processing, recently I read a paper on color image contrast enhancement, In that paper authors have implemented their method in HSV color model. In my work I have implemented that paper in MATLAB using RGB color model. I just want to know that whether it will make any difference or not?

HSV is one of the color space transforms that separate color from intensity. Depending on your application (my understanding is contrast enhancement is somehow related to image segmentation in your OP), there are several advantages of HSV:

Histogram equalization of a color image is suggested only on the intensity component;
The Saturation component is regarded as invariant in shadow (only intensity changes), so processing the image in HSV space is better suit for shadow removing. This post that tried to remove the shadow from an image is for your reference;
The Hue component makes the algorithm better immune and thus more robust to lighting variations. This feature makes it better fit in some medical applications such as the palm skin segmentation (under normal RGB space it is not easy to thresholding out the hand palm).
In color enhancement that you mentioned, if the colors in your image are not so strong, it would be quite difficult to detect colored objects. Adding some arbitrary value in the Image would not Help, since it would increase the brightness but not the contrast. After the conversion to HSV space you may increase the Saturation then convert it back to the RGB color space. The colors of interest may appear much better in the image which would facilitate the object recognition nicely. You can have a try by this means to adjust your image contrast.

You may also want to check this post for more discussions.


How to do simple extrapolation with Fourier transformation? I have 1024 sample points, and I would like to do really simple extrapolation using Fourier transformation. First I apply Fast fourier transformation on the data. My first intuition was that I just calculate the inverse fourier transformation on a larger interval. The result of the FFT is a collection of sinusoidal waves, so I expected that it's result is continuous at the end of the original interval.

I calculate the result with this formula:

$$
X[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{i 2 \pi n k}{N}}
$$

The extrapolation at the end of the interval looks like this:

I would like to avoid this "break" in the result.

The python source code is:

import numpy as np
import matplotlib.pyplot as plt

from scipy.fftpack import fft

# generate samples
data = np.zeros(1024)
for k in range(0, 1024):
    data[k] = np.sin(k/20 * np.pi*2)

# apply FFT on samples
waves = fft(data)

# calculate the new data points from the FFT result
res = np.zeros(2048)
for k in range(0, 2048):
    val = 0.0
    for n in range(0, len(waves)):
        # https://dsp.stackexchange.com/a/510/25943
        val += waves[n]*np.exp(1.j * 2*np.pi * n * k / len(waves)) / len(waves)
    res[k] = val.real

#plot the result
plt.plot(res)
plt.show()

Your computation is right. But your expectation is wrong.

The method that you use to extend the periodic sequence $\tilde{x}[n]$, of period $N$, outside of its base period  $x[n]$ for $0 \leq n <N$ into the range, say, $0 \leq n <2N$ will simply copy the base period $x[n]$ next to itself. Hence unless the last and first elements $x[N-1]$ and $x[0]$ are same (or very close) there will be an apparent mismatch (discontinuity), that you interpret as a glitch, at the boundary.

The base period $x[n]$ and  its $N$-point DFT sequence $X[k]$ are periodic by $N$ which means that $x[n+N] = x[n]$ and when extrapolated you should get the exact result. The inverse DFT that you implement is:

$$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} k n}$$

When evaluated for some $n = n_0 + N$ for the extrapolated samples $x[n_0 + N]$ where $n_0$ is inside the base period $0 \leq n < N$ you will have

$$x[n_0 + N] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} k (n_0 + N)} = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} kn_0} e^{j \frac{2 \pi}{N} N} = x[n_0] $$

as expected.


Visible light communication and LED modulation I work on a VLC system in which 16-PAM symbols with complex pilot symbols are as ifft input. I made the ifft output real and positive. Now how should i use these ofdm symbols to modulate LED? Their apmlitude or their power are used to be multiplied by time domain response time of LED?

OFDM in general is a complex baseband technology, which means it produces an I and a Q signal, which have to be mixed with a carrier wave and it's 90° shifted version, and added up before it can be used to excite a real-valued emitter (like an antenna or an LED driver).

By "I made the IFFT output real" you probably mean that you supplied the IFFT with a symmetric signal? That's nice, because it circumvents the need for a carrier frequency, but it also inherently halves your available bandwidth! Also, I can't fight the feeling that the discrete Fourier transform is not really what you're looking for if you just want to produce a real signal – the discrete cosine transform might be better suited for you.

The core question here is why you want to use OFDM in the first place.

The main motivation for OFDM comes from something that is highly inherent to the RF channel: Channels of high bandwidth aren't flat, so the equalizer for something like 20MHz of WiFi channel gets really ugly. Here, you can gain a lot by dividing your channel into orthogonal subchannels, and equalize them with a much simpler equalizer individually. I wasn't aware the same problem applies to the monochromatic VLC channel!


LTE Channel Estimation I am trying understand LTE channel estimation. I have multiple questions:

How much phase shift does a typical LTE channel have?
Can all that be corrected using pilots?
Where do we typically implement these algorithms? In a processor or FPGA? Because the OFDM system will be in an FPGA or similar platform.

FPGA is only user to verify the initial RTL which is specifically written to do each of the task listed below. Once the verification is done there is a tape out and everything is implemented in hardware accelerators as below. FPGA is wasteful of area and too power hungry which cannot be afforded in modem (low power device)

The LTE channel estimation is done as follows:

The FFT of the incoming samples at receiver is taken post cyclic
prefix removal. This is done in an FFT kernel in Hardware (specifically written RTL, not in  FPGA)
Then this result is fed to a memory where all the    subcarriers are
stored for this particular symbol (other symbols are    stored in
this memory as well).
Then firstly all positions at which    CRS    (Cell Reference Symbol)
is available are extracted for the symbol, the channel is estimated (MMSE or others) at the subcarriers corresponding to the CRS. This is done in a hardware accelerator. See figure below

Then there is an interpolation done in between the CRS subcarriers to
find the the channel estimate at non CRS subcarriers (wiener
filters), this is called interpolation in frequency. This is also
done in Hardware accelerators (i.e. RTL is written specifically for
this)
Now, CRS (per antenna port) is only transmitted in 2 symbols in a
slot(for most of the transmission modes). It  is transmitted in the
first and $5^{th}$ symbol of a slot for normal CP.
All the symbols for a full subframe  are buffered and frequency interpolation is done at all symbols in which CRS is transmitted as explained above (see figure below)

Then there is an interpolation done in time to get channel estimates
for all symbols in between the symbols containing CRS, to get a
channel estimate for all the subcarriers in the entire time frequency grid. This is also done in the Hardware accelerators
The CRS references positions are at every 6th subcarrier in frequency
(linked to the coherence bandwidth) and similarly in time linked to
coherence time

Once there are channel estimates available for all subcarriers  of 1 full subframe. It is sent for equalization/demodulation (this is also done in Hardware accelerators)

The reason for using Hardware accelerators over doing in Firmware is speed. The reason for doing in accelerators over FPGA is area and power. accelerators have a smaller footprint and more power efficient

To the first question, by itself there is no limitation that an LTE channel can have for the phase shift. But there are limitations as to what can be estimated unambiguously by the receiver. Consider this situation: The maximum frequency offset that an LTE system can estimate is +/-7.5KHz (because subcarrier spacing is 15 Khz) anything above this frequency offset will wrap around to within the above specified range. Now since frequency offset shows up as phase in the signal multiplied with channel estimate the ambiguity will remain in unwrapped phase.

As to how then we measure frequency offsets more than +/-7.5Khz, that's another story :)


What motion estimation techniques that could estimate the motion of every pixel? I'm trying to estimate a motion from a video, however based on my finding, those methods usually meant for large motion and sparse in nature.

Optical flow: Sparse and fail to detect sub-pixel motion

Phase-based: Only applicable to global shutter camera for accurate displacement

Welcome any suggestion/discussion you guys have in mind that would allow me to determine sub-pixel motion for every pixel that applicable to rolling shutter camera! Thanks in advance!

You can use Farneback to estimate dense flow.  In OpenCV, the function call is calcOpticalFlowFarneback().

Please note that the initial flow is important and you can use one of OPTFLOW_FARNEBACK_GAUSSIAN or OPTFLOW_USE_INITIAL_FLOW.

It depends on the depth of the pre calculation to estimate the segmentation before calculating the frame to frame flow.


Discrete-time Fourier transform I am a junior high school student who has a general fascination for electronics, programming, and the like. Recently, I have been learning about signal processing.

Unfortunately, I haven't done much calculus yet (forgive me), so I am a little fuzzy on things.

If you were to compute the DTFT of a signal, what would be the difference between a $\sin$ or $\cos$ representation of that signal?
With the DTFT I understand that the signal you input would be discrete in time, but how in the world can you achieve a continuous signal in the frequency domain?
This leads to my second question, which is: how is the DTFT useful? Where has it been used with most applications and why?

I would appreciate any help.

It's great that you are interested in signal processing at that early stage of your educational path.

The best path to get there is to read some introduction books on the topic. There are a lot of good and free online resources to get you started. [Note to the esteemed editor: good introduction books may a really good topic for a "sticky"].
I sometimes use

http://www.dsprelated.com/dspbooks
http://www.dspguide.com/pdfbook

One of the most important mathematical concepts you will need to get your arms around are “complex” numbers”. It’s clearly a misnomer since it’s really not that complicated and clearly makes nearly all engineering math a lot simpler. Another great free resource for all things math related is http://www.khanacademy.org and in this case specifically
http://www.khanacademy.org/video/complex-numbers--part-1?topic=core-algebra

Back to your first question: There are actually four different flavors of the Fourier Transform: Fourier Series (most likely to show up in high school), Fourier Transform, Discrete Fourier Transform and Discrete Fourier Series. All of them of them use a combination of both sine and cosine (or a complex exponential, which is essentially the same thing). You will need both.

Let’s say you calculate the sine and cosine Fourier coefficients of an input sine wave. (Under certain conditions) you’ll find that all Fourier coefficients will be zero except for one cosine and one sine coefficient. However, depending on the phase of the input sine wave you these two numbers will move around. You may get [0.707 0.707], or [1 0], or [0 -1], or [-0.866 0.5] etc. You will see that the sum of the squares of those two numbers will always be 1, but the actual values depend on the phase of the input sine wave.

If you want to deep dive, try this:
http://www.dsprelated.com/dspbooks/mdft/


correct for shift from butterworth filter? I have implemented a python butterworth filter but it seems to shift the data a lot, the x location of my data is very important to me for my analysis. Is there a way to correct for this shift precisely?

Does it shift all the data uniformly? Or does certain areas get shifted further/less far that others?

The delay or "shift" of an IIR butterworth filter varies with frequency.  If you want a constant delay that you can correct for in a simple manner, you need to instead use a linear phase FIR filter, not an IIR filter.


FFT of shifted image in Numpy I want to prove using numpy the theory of Fourier transforms in which translation in space corresponds to a shift in the phase domain (frequency domain remains constant).
I have generated these three 2d-numpy arrays showing translated versions of the same boxes:

I then use numpy to calculate the FFT of each of these images doing:
H1 = np.fft.fft2(im1)
H2 = np.fft.fft2(im2)
H3 = np.fft.fft2(im3)

Although the frequency and phase plots of the transformed arrays look very similar I perform a comparison (MSE) between them to see the differences.
This is the MSE plot for the frequency domain:

And this is the MSE plot for the phase domain:

I'm confused as I was expecting to see identical spectra in the frequency domain but different in the phase domain.
I would appreciate if someone can provide some insight to detect if the problem comes from a theoretical misunderstanding or a problem in the code.
Here is the code to replicate the results shown here.

In the code the OP is comparing the real and the imaginary components of the FFT in each to determine similarity. If a phase shift occurs then the real and imaginary components will also change. What the OP should be comparing is the magnitudes of the complex values in each sample, which should be invariant under such image shifts.
Confirming this; the OP had done the following to compare the FFT results:
plt.bar(['1-1','2-1','3-1'],[np.mean(np.square(H1.imag-H1.imag)),np.mean(np.square(H1.imag-H2.imag)),np.mean(np.square(H1.imag-H3.imag))])

resulting in:

If one does the following instead:
plt.bar(['1-1','2-1','3-1'],[np.mean(np.square(np.abs(H1)-np.abs(H1))),np.mean(np.square(np.abs(H1)-np.abs(H2))),np.mean(np.square(np.abs(H1)-np.abs(H3)))])

We get a result that we expect:

The equivalent magnitude is given by the translation property of the 2D DFT where:
$$x[k-m, l-n] \mapsto X[u,v]e^{-j2\pi (mu/M + nv/N)}$$
This is similar to the time shift property of the 1D DFT (which shows fundamentally what is occurring and can be extended to a 2D DFT):
Given 1D DFT as
$$x[n] \mapsto \sum_{n=0}^{N-1}x[n]W_N^{nk} = X[k]$$
Where $W_N^n=e^{-j2\pi n/N}$
From which:
$$x[n-m] \mapsto \sum_{n=0}^{N-1}x[n-m]W_N^{nk}$$
Substitute $u= n-m$
$$\sum_{u = -m}^{N-1-m}x[u]W_N^{(u+m)k}  = \sum_{u = -m}^{N-1-m}x[u]W_N^{uk}W_N^{mk}$$
$$=W_N^{mk}\sum_{u = -m}^{N-1-m}x[u]W_N^{uk}$$
for any k:
$$\sum_{u = -m}^{N-1-m}x[u]W_N^{uk} = \sum_{u=0}^{N-1}x[u]W_N^{uk}$$
Therefore
$$x[n-m] \mapsto W_N^{mk}X[k]$$
And $W_N^{mk}$ as $e^{-j2\pi mk/N}$ only modifies phase, not magnitude
Also the DFT is circularly periodic in the time domain, so any such rotation must done in circular fashion;  x[mod(n-m,M)] = x [n]


Attenuating peaks of a signal in a specific direction Consider a segment of a time-domain signal shown in the following graph and the two marked peaks, the first one facing downwards (in red) and the second one facing upwards (in black).

I would like to attenuate the peak facing downwards only. Unfortunately, both peaks seem to have the same frequency, so a linear filter that attenuates that frequency will attenuate both peaks.
Is there a kind of a filter (perhaps a non-linear filter) that can attenuate peaks that face downward and keep peaks that faces upward (or vice-versa)?

My immediate approach to this: try something that I'd call a running-minimum filter:

Look at a fixed amount of $N$ samples around the current sample. If the current sample is significantly below that, modify the sample to be closer to a low-pass variant of the original signal.


MATLAB's $\tt filter()$ function for complex valued data I want to use the filter() function to implement an FIR filter having complex data values as the input to it. For real, I used this function as filter(A,1,B).

Question :

Would there be a difference for the corresponding implementation when A and B are in complex domain?

As @msm said, there's no difference. You can refer to Matlab's documentation here: https://www.mathworks.com/help/matlab/ref/filter.html#inputarg_x which says "Complex Number Support: Yes".


What's this discrete latch operation called? I have a picture in my head of a specific DSP operation, but I don't know what it's called or how to google for it. I'm pretty sure it's a fairly common operation, however I have no idea what to call it.

Consider a discrete event that occurs at time interval (Excuse the ascii art):

                          |
              |           |
x       |     |  |        |
     _  |  _  |  |  _  _  |  _  _
     0  1  2  3  4  5  6  7  8  9
                     t

What I have is two vectors descibing this waveform.

x = [2,3,2,4]
t = [1,3,4,7]

What I want to do, is 'latch' the most recent signal over time, like this:

                          |  |  |
              |           |  |  |
x       |  |  |  |  |  |  |  |  |
     _  |  |  |  |  |  |  |  |  |
     0  1  2  3  4  5  6  7  8  9
                     t

Which would produce this waveform:

x = [0,2,2,3,2,2,2,4,4,4]
t = [0,1,2,3,4,5,6,7,8,9]

In electrical terms, this would be a memory element which remembers the most recent signal received. It also would be some crude interpolate function, however all the interpolates I know about (linear, cubic) fill in the in-between with more complicated functions.

I'd like to recreate this behaviour in Matlab but don't know how to google for what I want. So what's this called?

I figured this out, it's called a Sample and Hold. The matlab page on it talks about how to implement one in matlab.

In summary, an interpolate function is used to create it. There are a couple of interpolate options, like "next", "previous" and "nearest". The "extrap" is needed to extrapolate the waveform for times 8 and 9. Continuing my example:

x     = [2,3,2,4]
t     = [1,3,4,7]
t_new = 1:9

x_new = interp1(t,x,t_new,"previous","extrap")

Produces

x_new =
     2   2   3   2   2   2   4   4   4


Why is PSNR used for image quality metrics instead of SNR? I have been trying to determine why is PSNR used instead of SNR for image quality metrics?

I have read through several books and most just say that PSNR is used instead of SNR, but don't give a reason as to why.

I also read these posts:

Signal quality metric

and

Difference between SNR and PSNR

But I still feel like they don't explain clearly why PSNR is used instead of SNR.

I know that PSNR measures how the peak signal intensity compares to the mean squared error whereas SNR measures signal power compared to noise, but this doesn't help me reach an intuitive conclusion as to why to use PSNR instead of SNR.

Can someone please explain this?

Consider 8-bit quantization. This gives a mean squared error of about 0.577 least significant bit (LSB). Quantizing two images of different average intensities gives two different SNR, because that metric is normalized against the average intensity of the input image, but identical PSNR, because it is normalized against the maximum pixel value of 255. When describing the quality of something like a compression method, it is better to use a metric that is (to an extent) independent of the source material, because it carries meaning even when communicated without providing the input image.


List of Fourier coefficientes to simulate musical instruments I'm trying to find a list of Fourier Coefficients (frequencies/amplitudes, ie, spectral data) to simulate different musical instruments using additive synthesis. Are these data published somewhere? This is for a Math class: I just would like to add more frequencies and listen how the sound generated compares with the real sound.

Thanks!

look up Additive Synthesis somewhere.  those Fourier coefficients are not constants.  the musical note is not normally a perfectly periodic function in time.  it is what i like to call "quasi-periodic".  so this is not exactly true (given some period $P$):

$$ x(t) = x(t+P) \quad \text{for all } t $$

but in the neighborhood of some given time $t_0$ this is approximately true

$$ x(t) \approx x(t+P(t_0)) \quad \text{for } t \approx t_0 $$

so the exact Fourier series for $x(t)$

$$ x(t) = \sum\limits_{n=-\infty}^{+\infty} c_n \ e^{j 2 \pi n f_0 t} \quad f_0 \triangleq \frac{1}{P} $$

is not precisely true, but more like

$$ x(t) = \sum\limits_{n=-\infty}^{+\infty} c_n(t) \ e^{j \theta_n(t)} $$

where

$$ \theta_n(t) = \int\limits_{0}^{t} 2 \pi f_n(u) du + \theta_n(0) \quad f_n(t) \triangleq \frac{n}{P(t)} $$

so what you're looking for is a collection of envelopes, rather than a collection of constant coefficients, for some musical instrument.

in a previous life, living on an older computer (that might not power up now) and maybe even on some 3.5 inch floppy disks, i have some of that additive synth data.  actually, i am not sure where to get it all digested.  within the context of Wavetable Synthesis, i know how i would go about extracting it from "samples" of notes (like what you would get from a sample library like from Synthogy ).


Looking for a Mathematical Derivation for the Energy Formula in Continuous-time Domain I have just started my signal and system course and I would like to know how we derive the corresponding formula for the energy of a continuous-time signal $x(t)$ over an interval $[t_{1},t_{2}]$ :

$$E_{x}=\displaystyle\int_{t_{1}}^{t_{2}}|x(t)|^{2}\;\text{d}t$$

That is to say why is the signal energy the area under the region of the square of $x(t)$. For instance, I can think of $x(t):=t$ with $t$ being geometrically the first bisector and so by squaring it we obtain a parabola so this would mean the energy is the area under this parabola but why is this true?

The formula for the energy of a signal over an interval isn't a derivation.  It's a definition.  You do that integration and you call the result the signal's energy.
The reason that in the mathematical world we make that definition is because (A) it makes the math handy, and (B) for a preponderance (but by no means all) of the means of signal transmission in the real world, the actual, physical, energy of a signal over an interval really is proportional to its amplitude squared.
So if you launch a radio signal into free space with an antenna, you couple into free space where the actual physical energy of the signal is proportional to its amplitude squared.  When you launch a sonar signal into the water with a transducer, the actual physical energy of the signal is proportional to its amplitude squared.  When you launch an electronic signal into a transmission line, the actual energy of the forward-traveling portion of the resulting wave is proportional to the amplitude squared.  Etc.
So -- just accept it because it makes the math handy, and be happy that the real world happens to fit with an expression that also makes the math handy.


What is the link between equation of a continuous signal versus equation of its sampled form? I have just started DSP self-learning. I am a little confused by this 'end of the chapter' exercise question from Chapter 2 of "Understanding Digital signal processing" 3rd edition by Richard G. Lyons.

  Consider a continuous time-domain sine wave defined by

  $$x(t)=\cos(4000\pi t)$$

  that was sampled to produce the discrete sine wave sequence defined by

  $$x(n)=\cos(n \pi/2)$$

  What is the sample rate ($f_s$ measured in $\textrm{Hz}$) that would result in the sequence $x(n)$?

There are a few things that I have not understood about this question:

Why is the argument of the cosine function so different after sampling?
I am not sure why the argument of the cosine after sampling does not include $t$ so that we can get $nT$ where $T$ is the sample time.
What determined the argument of the cosine after sampling?
Can I determine the frequency of the sampled signal from the argument of the cosine in the $x(n)$ equation?
And of course a cheeky one, the answer to the book question please :)

The sequence $x[n]$ equals the continuous-time signal $x(t)$ sampled at $t=nT$, where $T=1/f_s$ is the sampling period:

$$x[n]=x(nT)=\cos(4000\pi nT)=\cos(4000\pi n/f_s)\tag{1}$$

So if $x[n]=\cos(n\pi/2)$ you just have to compare this argument to the argument in $(1)$ to figure out what $f_s$ is.

The normalized frequency of the discrete-time signal $x[n]$ is $\omega_0=\pi/2$ (in radians). It is related to the frequency $f$ in Hertz by

$$\omega_0=\frac{2\pi f}{f_s}$$


I/Q Signal Power Calculation I am trying to understand the basics of IQ modulation. The amplitude of a demodulated IQ signal is given as $\sqrt{I^2 + Q^2}$ which is easy to understand from the phasor diagram.

Yet, I don't understand the IQ power: ${I^2 + Q^2}$ . Where is the load impedance in this expression ? I'd expect something like

$Z_{load} = 50 \Omega $

$P = \frac{I^2 + Q^2}{Z_{load}}$

What is missing in my understanding ?

Nothing is missing in your understanding. When we calculate power that way, it is assuming a $1 \Omega$ load and that amplitudes are given in volt, but almost nobody ever says that even in textbooks. Unless it is explicitly stated otherwise, I think it is safe to assume this.


Do DSP principles apply in the same way to both Audio and EM spectrums? All  of my DSP experience is very much centred in the EM spectrum and I've only ever processed EM-type signals in the past.

I'm now doing quite a bit of work in ultrasound and the audio spectrum, but a lot of the techniques I'm using are used in the EM mainly. I'm reading a lot of EM papers about DSP at much higher frequencies. Can I expect that the principles behind these techniques will behave in the same way, even though the medium of transmission is different?

Many of the ideas will be similar, but some assumptions may no longer apply. For instance many signals in EM can be treated as Narrowband (narrow bandwidth with respect to the centre frequency), in Ultrasound you may run into Broadband signals more easily. Another problem you may run into is the speed of propagation varies as a function of frequency.

Exactly what difference you encounter will depend on your experiences with EM and what exactly your are doing in ultrasound.


Signal that destabilises a positive feedback path with nonunity feedback I am just starting to learn about feedback systems. There are a few questions that have me stumped, this being one of them. Could someone please help? Thanks!
Question: Find a signal that destabilizes a system with a delayed positive feedback path and nonunity feedback gain.
My understanding so far is that we are reintroducing the output as an input into the system with a positive coefficient. So if the input is positive, the signal with be amplified. Nonunity feedback gain means only a fraction of the output is being fed back into the system. Is that correct so far?
Now we want to find a signal that is destabilised by this feedback system. That means the original signal should be stable.
A sinusoidal wouldn't work, because the negative inputs cancel out the positives. The square of a sinusoidal is already unstable. How about the absolute value of a sinusoidal? Would that work?

I am going to provide some initial information as simply as I can that should be helpful in the OP's attempt to get a basic understanding of control systems. Further study in Laplace transforms (and z-transforms for discrete-time systems) will be helpful in these frequency domain modelling approaches to control systems.
First understand the concept of a transfer function as the frequency response for a basic two port (one input / one output) linear system; the transfer function defines the magnitude and phase versus frequency of the output compared to the input.  This is often given in terms of the variable "s", but to avoid going down that path I am going to present it as a function of frequency only ($H(f)$).  Below shows an example system with a plot of the magnitude and phase versus frequency as given by the $H(f)$ equation. This shows us for any particular frequency at the input, what the amplitude of it will be at the output and what the phase shift will be between the input and output: every frequency will be transformed by the system by changing its amplitude and phase shift only, as long as it is a linear system no other frequencies will be created so this is a nice compact description that we will ultimately use to determine stability in a control system.

Now consider the simplest case of a system with feedback as diagrammed below with input given as $X$ and output given as $Y$ with the output fed back and subtracted from the input:

The relationship between $Y$ and $X$ is easy to show with simple algebra, with our goal to describe a new transfer function of this new closed-loop system. The new transfer function would be given as $Y/X$ (output divided by input). I don't show it below to not clutter the simplicity of the equation but $X$ and $Y$ are also functions of frequency as $X(f)$ and $Y(f)$, so would be the Fourier Transforms of time domain signals:
$$(X-Y)H(f) = Y$$
$$\frac{Y}{X} = \frac{H(f)}{1+H(f)}$$
The above equation would adequately describe our new closed-loop system as long as it follow the Nyquist Stability Criterion. To really show this completely to cover all cases, the OP needs to first understand the Laplace Transform and the generalized transfer function as a function of s (as in $H(s)$), but here in this simple introduction we can reduce it to a condition that is covered in most cases: that is if the gain of the system given by $|H(f)|$ is less than one before the phase of the system given by $\angle (H(f))$ exceeds 180°, the system will be stable and will have its own closed loop transfer function given by the expression above. We can see this through an iteration of a simple case with H(f) = -.8 at a particular frequency (gain = 0.8 and phase = 180°).
Consider the feedback system as an iteration with the input starting at 1: the output will be -.8, which fed back to the input as a subtraction will add 0.8 to the input, resulting in 1.8 at the output. This is then multiplied by -.8 resulting in 1.44 being added to the input resulting in 2.44 at the output...
With closer inspection we see that this is a geometric series:
$$1 + \alpha + \alpha^2 + \ldots = \sum_{n=0}^\infty(\alpha^n) = \frac{1}{1-\alpha}$$
Which converges only if $-1 < \alpha < 1$.
Where here $\alpha$ was 0.8 (due to the negative feedback in the loop) which would converge to $1/(1-.8) = 5.
If the gain of $H(f)$ exceeded one with the phase 180°, we would see how the same geometric series would not converge, as an example of an unstable loop.
This gives a very basic simplified explanation to help aid in an intuitive understanding, but further proper study would be to understand the Laplace Transform and complex systems, and the meaning of poles and zeros in the s-plane, and how a stable system must have all of its poles in the left-half plane, and ultimately the full details of the Nyquist Stability Criterion.


Why isn't the Hartley Transform more widely used? I was doing a bit of research on FFT and I came across this paper. According to it, computing the Fourier Transform using FHT is significantly faster than other methods, even for complex data. Yet, from what I gathered so far, the FHT method is rarely used. Why is that?

Basically, it isn't true.  A Fourier transform making use of the symmetries of real-valued original data is just as computationally efficient.  Also you can "pack" two independent sets of real data into one complex transform (one in the real part, one in the imaginary), then use symmetry relations to recover the respective two transforms.

The Hartley transform butterfly operation requires working with more input and output operands than the Fourier transform butterfly once the "twiddle factor" becomes smaller.  Basically you have have two forward indexed data sets and two backward indexed data sets flowing into two butterflies in order to get two forward and two backward indexed data sets out in-place.

So all in all, the number of operations does not become less as soon as your algorithms make use of the "real input data" criterion, and the indexing logic is a lot less transparent.  It's similar with the convolution theorem: applying it to Hartley transform data does not really save operations, you need to work with forward/backward address pairs and the algebraic operations are more tenuously related to the analytic relations of the Fourier transform than if you work with FFT right away (possibly including the Hermitian properties of a real-valued original transform pair).

The data representation is nicer if you want to work with only the non-redundant elements of a transform of purely real-valued original data (packing the real-valued $f_0$ and $f_{N/2}$ into one complex data element is ugly when storing non-redundant FFT data).  Indeed, considering the Hartley transform as a dual storage representation of FFT (basically even for complex original data, $H(k) = \frac12 (F(k)+F^*(-k))$ if I remember correctly) at each stage of computation is one of the more straightforward ways to figure out a DHT step by step.


What Resources Are Recommended for an Introduction to Signal Processing (DSP)? What are some recommended resources (books, tutorials, lectures, etc.) on digital signal processing, and how to begin working with it on a technical level?

My recommendation in terms of text books is Rick Lyons's Understanding DSP.  My review of the latest edition is here.

I, and many others from the ${\tt comp.dsp}$ community and elsewhere, have helped Rick revise parts of the text since the first edition.

For self-study, I know of no better book.

As an on-line, free resource, I recommend Steve Smith's book.  Personally, I prefer Rick's style, but Steve's book as the advantage of online accessibility (and the online version is free!).

Edit:

Rick sent me some feedback that I thought I'd share here:

For your colleagues that have a copy of my DSP book,
I'll be happy to send them the errata for my book.  All they have
to do is send me an E-mail telling me (1) The Edition Number,
and (2) the Printing Number of their copy of the book.  The Printing
Number can be found on the page just before the 'Dedication' page.
My E-mail address is:
R.Lyons [at] ieee.org

I recommend that your colleagues have a look at:
http://www.redcedar.com/learndsp.htm

Rick also gave me a long list of online DSP references.  There are way too many to put here. I will see about setting up a GoogleDocs version and re-post here later.


Wavelet central frequency What is the definition of the central frequency of a Wavelet and the bandwidth of a Wavelet? How do they affect the resultant waveform of the wavelet?

Thanks.

It's not formally defined, however I've seen $\int_{-\infty}^\infty \omega |\Psi(\omega)|^2 d\omega$ be used in the literature (i.e. the centre of mass of $|\Psi(\omega)|^2$), where $\Psi$ is the Fourier transform of the wavelet $\psi$.