Is there a condition for bandpass sampling? Consider a signal that has lowest frequency component $F_l$ and highest frequency component $F_h$.
According to the theory of bandpass sampling, this signal can be sampled and succesfully recoverd if sampled at a frequency ($F_s$) twice the difference between highest and lowest frequency.

that is, $F_s=2 \cdot \left(F_h-F_l\right)$.

My doubt is , Is there a condition or restriction on $F_h$ and $F_l$ for bandpass sampling?

It's generally not true that a band pass signal can be sampled and recovered without error if $f_s>2B$ is satisfied, where $B=f_h-f_l$ is the signal's bandwidth. This condition is just necessary but not sufficient.

You have to make sure that the aliased spectra do not overlap. This results in the following condition on the sampling frequency:

$$\frac{2f_h}{n+1}<f_s<\frac{2f_l}{n}\tag{1},\quad n=0,1,\ldots$$

If Eq. $(1)$ can only be satisfied for $n=0$, you get the familiar Nyquist sampling theorem where $f_s$ must be greater than twice the highest frequency of the signal: $f_s>2f_h$. The lowest possible sampling frequency is obtained for the largest integer $n$ such that $(1)$ is still satisfied. This maximum value of $n$ is given by

$$n_{max}=\left\lfloor\frac{f_l}{B}\right\rfloor\tag{2}$$

where $B=f_h-f_l$ is the bandwidth. Eq. $(2)$ is obtained from $(1)$ by computing the largest value of $n$ such that the left-most term is still less than the right-most term.

As an example, take a band pass signal with $f_l=10\text{ kHz}$ and $f_h=30\text{ kHz}$. From $(2)$, Eq. $(1)$ can only be satisfied for $n=0$, so you have to choose $f_s>2f_h=60\text{ kHz}$, which is a stricter condition than just requiring $f_s>2B=40\text{ kHz}$. However, if $f_l=10\text{ kHz}$ and $f_h=14\text{ kHz}$, the maximum value of $n$ for which $(1)$ is satisfied is $n=2$, which gives you a range $2\cdot 14/3=9.33<f_s<10$ of possible sampling frequencies, all of which of course satisfy $f_s>2B=8\text{ kHz}$.


Servo motor analysis I'm studying a mathematical behaviour of a servo motor and I need some help to understand it.

The output signal is $\beta(t)$, representing the angle rotated by the axis at instant t, in relation to the equilibrium position.
On the servomotor there are two torques:

the first torque is produced by the electric current $i(t)$ that does through the motor. It is proportional to the electric current being given by $Ki(t)$.
the other torque is given by friction and is given by $-2 \beta'(t)$ where $\beta'(t)$ is the derivative of $\beta(t)$, meaning it represents the angular velocity.

The feedback control system is given by:

Where Motor represents the differential equation
$ 4\beta''(t)=100i(t) - 2 \beta'(t) $

And $\alpha(t)$ is the input angle that give the desired position (angle) for the servomotor.

$B$ is unknown but we will assume $A = -1$ so that if the input is constant we'll have the output equal to the input (unitary static gain).

The transfer function of the system is:

$$H(s)= \frac{100}{4s^2+(2-100B)s-100A}$$

By this expression we conclude that our system is a second order system with no zeros.
By controlling the value of $B$ we will control where the poles are in our system and by that the type of system we get. More specifically:

If $B<-0.38$ we have two real poles (and an overdamped system)
If $B=-0.38$ we have o double pole (critically damped system)
If $-0.38<B<0.02$ we have two complex conjugated poles (underdamped system).

So until this I don't have questions and everything makes sense to me.
The next affirmation is what leaves me doubtful:

By controlling the value of $B$ we are controlling the friction of the system.

No more explanation is given to this and I'm kinda worried if my interpretation of this fact is correct or not. Here it goes:

We have a 2nd order system without zeros. By controlling $B$ we will control the type of poles our system will have. For low values (in module) of $B$ we have an underdamped system. For $B=-0.38$ the system is critically damped. For $B<-0.38$ (meaning, large absolute values of $B$) the system will be overdamped. The system damping is related to the existence of friction in the motor. For low absolute values we will have little friction and so the system will oscillate a lot until it reaches its final value. As we increase $B$, we will increase friction which will oppose to the natural oscillatory tendency of the system, making it monotonous until it reaches its final values. We have to be careful however because as we increase the friction the system will also take more time to reach the value (it has to "beat" the friction).

Is this interpretation correct? I'm not sure if this is the correct relationship between friction and the type of system I get? I'm also curious if there is any more implication or consequence of this relationship. Can anyone clarify me please?

Thanks!

Yes your interpretation is correct. Your model of the servo motor creates a second order differential equation:
$$I \beta'' + b \beta' + k \beta = \tau_{ext}$$

where $\beta$ is the position (in radian angle) , $\tau$ is the externally applied torque (given by the input current and the mechanical load if it exists and is being modeled), $I$ is the moment of inertia of the rotating part, $b$ represents the effect of frictive forces which depend on angular velocity such as your damping friction effect, and $k$ is the effect of position dependent torsional-spring force if it exists. In your system you only consider the angular velocity dependent frictive force and ignore the (non-existing indeed) spring force hence $k=0$. This yields the differential equation

$$I \beta'' + b \beta' = \tau_{ext}$$

The system can also be analysed by Laplace transforming the LCCDE, at this point I'm not sure how you reached your $H(s)$, but assuming it's correct then it can be seen that the coefficient $2-100B$ models the frictive force strength and as you have stated; by adjusting $B$ you will adjust the type of solutions from oscillatory underdamped to monotonic overdamped.

And the physical explanation of the result is also in line with your intuition; with more friction the motor will hardly oscillate but will be slow and inefficient while with less friction the motor will respond faster but will overshoot and osccillate. The critical damping is the optimal solution which reaches the equilibrium as fast fast possible without oscillations. That could be relaxed a little for practical reasons and you can allow some tiny oscillations to improve the dynamic performance considerably, of course provided that a minute amount of overshoot is tolareable.


How do you estimate necessary dynamic range? Say I have a signal (1D in my case, but could be an image, etc.) that has been captured adequately on nice equipment with sufficient resolution and range to get all the details of the signal that are important. I want to estimate what the minimum dynamic range to capture the relevant details of the signal would be.

In particular I want to compare the dynamic range of two signals.

My idea was to take the top percentile of the absolute value of the stored signal, divided by the bottom percentile. In python, for example:

DRestimate = np.percentile(y,99)/np.percentile(y,1)

Obviously this depends on subjectively picking a good percentile for your application, but it does allow me to get a comparison between two signals. So, the question is, is there a standard, or better way to do this?

Well, as you guessed, it depends!

For example, assume you're observing a multi-carrier communication signal. It has probability of 0.1% of symbol times to land in upper the percentile you "cut off".

Let's say that in that case, you get a bit error rate of 25% whenever signal is cut off.

Now, you get a BER –due to cutting off– of 2.5·10⁻³. That's pretty much for some systems, and nearly nothing for others.

So, define this from what your measurement describes, and generate a curve of "error you introduce" over "cutting off at a specific level"; let's call that function $e(l)$ ("error of level").

Then, calculate how often that level $l$ or anything lower occurs; that's a cumulative distribution function of the signal amplitude! It can be easily estimated from an observation. Because it's a probability distribution, we'll call it $F_L(l)$ of the random variable $L$, the observed level.

$\tilde F_L(l) = 1 - F_L(l)$ is just the "counterpart": How often does a value above $l$ occur?

Now, you can find the cumulative error of a cutoff-level by just defining a function

$$C(l) = \tilde F_L(l) e(l)$$

and evaluate it for some different $l$. (Do a plot of this, really.)

Then, calculate the percentage of broken measurements you're willing to accept, and compare your $C(l)$ to it.


SAR spotlight raw signal database For practicing spotlight processing, i am looking for I/Q baseband signal from a SAR airborne or spaceborne system. In my research I have only found processed image from :

eoweb geoportal https://eoweb.dlr.de/egp/
Airbus Terrasar X https://www.intelligence-airbusds.com/imagery/constellation/

Could anyone give me a link for raw I/Q signal (and associated metadata)?
Thanks a lot

This website offers phase history (match filtered) data: https://www.sdms.afrl.af.mil/index.php?collection=gotcha. There are other datasets on the same website, some of which may have raw I/Q.


Relative merits of fixed point vs floating point computation? I have a digital signal processing system that operates on a fast x86 machine using double precision floating point numbers.
It occurred to me that I am not really using the huge dynamic range of the floating point representation -- all quantities fit easily in the range ±32768.

My question: is it possible that switching to fixed point computations would provide a benefit in numerical accuracy (high priority) or computation time (low priority)?

Of course, the answer depends on how many bits are available for a fixed-point calculation.  How many bits of precision do typical fixed-point systems utilize?  Is it possible to efficiently do fixed-point calculations, with, say, 64 bits (16 bits integer part, 48 bits fractional part) on x86-64?

I had always thought that fixed-point calculations were used only in situations where CPU power is limited--does it make sense to use fixed-point calculations when CPU power is no concern?

The numerical accuracy of integers will only be better than the numerical accuracy of floats if the integer resolution is better.  Doubles have 52 fractional bits, so double-precision floats have a resolution worse than integers at around $2^{52}$, which is much larger than 32768 ($2^{15}$).  So, no, the numerical accuracy will not be better if you go to integers.

The second issue is speed.  The answer is: it depends on the hardware.  If you are running your program on a digital signal processor that has multiple fixed point multiply/accumulate cores then yes, it will be far faster in fixed point.  On an x86 chip, on the other hand, it will actually probably be slower in fixed point.  I did exactly what you are talking about once and saw my execution times increase.

After doing some searches on the internet I found out that that is common.  The reason is because it has a dedicated floating-point processor that isn't doing anything when you transition to fixed point, while the fixed-point hardware is shared with the regular fixed-point action, such as pointer arithmetic.

If you want to speed up the processing, the way to do it is to change from double-precision floats to single-precision floats.  That should produce a significant increase in speed.  That would, of course, reduce your numerical accuracy though.


Pedestrian counting algorithm Currently I am developing a pedestrian counter project (using OpenCV+QT on Linux). My idea about the approach is:

Capture Frames
Do Background Subtraction
clear noises (erode, dilate)
find blobs (cvBlobslib) - foreground objects
For each blob, set ROI and search for pedestrians (LBP with detectMultiScale) in these blobs (for better performance)
For each found pedestrian do a nested upper body search(Not sure) (better reliability)
If same pedestrian is found on continuing frames (3-4 frames maybe) - add that area to camshift and track - mark as pedestrian
Exclude camshift tracked areas from blob detection for next frames
If a pedestrian crosses a line increment number

I want to check if I am on the right track. Do you have any suggestions on how to improve my approach? If somebody worked on something similar, I would appreciate any useful tips, resources (and criticisms) on this problem.

I think what you are asking is about the feasibility of your pedestrian algorithm.

There are two general strategies for this kind of problems:

(bottom-to-top) Consider it as a pure detection problem, where in
each frame you detect only pedestrians. Once you detect them, a)
counting their number in a frame is fairly easy; and b) tracking any
of them in consecutive frames is also easy. Therefore, you solve everything.
(top-to-bottom) Consider it as a action recognition problem, where you detect whether a ROI in consecutive frames is a pedestrian or not according its actions. Once you solve this problem, then you solve the detection and tracking problems simultaneously.

Your heuristic algorithm is in the first category. I donot want to discourage you, but you might miss the key point that how to detect a pedestrian. Because real data can be more complicated than what you thought. For example, if this is a subway camera, whose frames maybe full of pedestrians and thus removing background or detecting blobs does not help at all. In this case, it maybe more reasonable to use face detection and face recognition algorithms to solve the problem, because if you find a face you find a pedestrian. On the other hand, depending on your definition of a pedestrian, it maybe true that not everyone appeared in a frame should be treated as a pedestrian. In this case, it might be reasonable to use action recognition algorithms (the second category), where you can explicitly define a pedestrian based on his or her behaviors.

Here are some tips based on my experience:

Stick on what you know and what you can easily pick up. Donot invest
your time on something is fancy but requiring many backgrounds you
donot have. Trust me, all these algorithms will be good in some case
but bad in some others. Therefore, the first thing is to make
something work, no matter how good or bad it is.
Know more about your data and then determine your method. A general description about a problem is insufficient in many cases.
If you want to demonstrate your idea,  it is better to use MATLAB and build a prototype.


Additive Gaussian white noise - bandwidth The average power at output when additive Gaussian white noise is present is given by the formula below:

$$P_0 = \int\limits_{-B}^{B} \frac N2\ df ,\quad \text{with  $B$ being the bandwidth.}$$

My question is when the bandwidth is in terms of $\textrm{rad/s}$, what would be the formula?
The one stated divided by $2\pi$ or another formula?

Just to be clear about the units:

\begin{align}
\Omega &= 2\pi F, \quad \text{with $\Omega$ in } \rm\underline{radians/sec}\text{ and $F$ in } \rm\underline{cycles/sec} \text{ or }\rm\underline{Hz}\tag{1}\\
\omega &= 2\pi f, \quad \text{with $\omega$ in } \rm\underline{radians/samples}\text{ and $f$ in } \rm\underline{cycle/sample}\tag{2}\\
\omega&=\Omega T_s\quad \text{and}\quad f=F/F_s\tag{3}
\end{align}
Read @Dilip Sarwate's comments; what should happen there is a change of variables.
\begin{align}
\omega &= 2\pi f\Longrightarrow d\omega = 2\pi df\iff df=\frac{d\omega}{2\pi}\tag{4}\\
f&=\pm B \Longrightarrow \omega= \pm 2\pi B\tag{5}\\
\end{align}
Using $(4)$ and $(5)$ you then have the following results:
$$
P_0=\int_{-B}^{B} \frac N2 df = NB \Longrightarrow P_0= \int_{-2\pi B}^{2\pi B} \frac N2 \cdot \frac{1}{2\pi}d\omega=\frac{1}{2\pi}\int_{-2\pi B}^{2\pi B} \frac N2 d\omega = NB\tag{6}\\
$$


Confused on Pulse Filter Bandwidth and Symbol Rate Relation ***Still need help in 2021 - not fully clear still on Jan 20th ***
I am confused about the relation between Sinc and Rectangle transform pair and how that relates the Bandwidth of Pulses, Bandwidth of Zero-ISI Filter and the Symbol Rate.
Using :

$T_s$ = Symbol Time Period
$R_s$ = Symbol Rate
$Sinc(t/T_s)$

Logic of Bandwidth of Sinc or Raised Cosine roll off = 0

Time-Domain Sinc : Every lobe is width $T_s$ but main lobe is $2T_s$ because its symmetrical on zero
Frequency-Domain Rectangle : Positive frequency Width of $\frac{1}{T_s}$ (But its Symmetrical on zero and so BW include negative frequency)
So...BW $= 2\frac{1}{T_s} =2 R_s$

Logic of Bandwidth Rectangle

Time-Domain Rectangle : Width $T_s$ (symmetrical on zero)
Frequency-Domain Sinc : Every lobe width is $\frac{1}{T_s}$ but main lobe is $\frac{2}{T_s}$ because its symmetrical on zero
So... BW $= \frac{2}{T_s} = 2R_s$

But I have read that's not right, the Zero ISI Filter BW = $R_s/2$ when the roll off factor  is zero. Why has my logic above got this wrong?
***Still need help in 2021 - not fully clear still on Jan 20th ***

Nyquist says that you can send up to twice as many pulses (symbols) per second as the channel bandwidth $B$ with zero ISI, so you need $R_s \leq 2B$. That is all there is to it. The sinc pulse has zero excess bandwidth so the bandwidth of the signal is equal to the symbol rate, $B_s = R_s$.
About sampling, Nyquist says that you need to sample at least twice as fast as the bandwidth of the signal, so you need $F_s \geq 2B_s$, where $F_s$ is the sampling rate.
Edit
To more directly address the points in the question...

Bandwidth of pulses: the FT of a perfect rectangular pulse is a sinc function which has infinite bandwidth.

Bandwidth of zero-ISI filter: the FT of a sinc function is a perfect rectangular pulse. If the zeros in the sinc function are at $kT_s$, where $k \in (-\infty, \infty)$, then the bandwidth is $\frac{1}{T_s}=R_s$. In case the time domain signal is purely real, then the FT will be symmetric meaning that the negative frequencies give you no more information then the positive frequencies. This cuts the bandwidth in half to $\frac{R_s}{2}$.


DSP with 8+ inputs for simple beamforming Does anyone know of a DSP that is cheap and has 8+ inputs (8 I2s mics) and that can perform beam forming calculations? I'm super new to circuit design so preferably one thats beginner friendly.

Thanks

edit: We are students who were kind of thrown onto this project for a class with little prior experience and one of the requirements is more than 8 microphones.

Details on our actual school project:
This is a sponsored senior design project, where a professor and one of his PhD students asked us to build a large (8+) microphone array that can be worn on a person's body. This array would collect data and store it in some data storage device that our sponsors can later access and run algorithms on. How we accomplish this was left open. The larger the array the better, our budget is unknown but likely $100-200 range.

For our class we have to add additional requirements so we decided to attempt to use the data collected and run a beamforming algorithm on it to amplify noise in a particular direction relative to the person wearing the device. Essentially our demo would be someone wearing the device and pointing it at different sound sources and hearing the amplified audio through headphones.

Essentially for the class we need to demo the device well but also build the device in such a way that our sponsors could use it for their own purposes after the class is over.

We all have some 300 level DSP knowledge but it is all "textbook" knowledge and no hands-on hardware knowledge so we are attempting to learn it as we go.

So far we have a microphone picked out here (ofc this can be changed) and are attempting to determine what DSP to use so that we can play additional parts of the device such as how we are going to transfer data to the headphones and data storage. It is also important to have some hardware component in the project which is why Arduinos and Raspberry Pi's are discouraged generally.

8+ high quality inputs for beamfroming and "...beginner friendly..." are competing requirements.
There will soon be the audioinjector which attaches on Raspberry Pi and would be a relatively pain-free option for what you mention. You can program the Raspberry in a number of different ways including "low-level" C to high level Python or other.
Another option would be to look for a multi-channel USB sound card which is Linux-friendly and attach it, again, to the Raspberry Pi (or small form factor pc). An older iteration of this one was just great for this purpose. And you get low noise amps, phantom power for mics and an internal DSP engine for basic preprocessing (although that was not accessible to the linux driver, at least a few years ago). This option again would give you the ability to program the multi-channel aspect via established hardware and software interfaces. Also, you might want to have a look at Jack and Portaudio.
Hope this helps.
EDIT:
I appreciate the additional information. As it is in narrative form, I am going to try to unpick the requirements. You say (emphasis is mine):

...a professor and one of his PhD students asked us to build a large (8+) microphone array that can be worn on a person's body. This array would collect data and store it in some data storage device that our sponsors can later access and run algorithms on. How we accomplish this was left open. The larger the array the better, our budget is unknown but likely $100-200 range.
For our class we have to add additional requirements so we decided to attempt to use the data collected and run a beamforming algorithm on it to amplify noise in a particular direction relative to the person wearing the device. Essentially our demo would be someone wearing the device and pointing it at different sound sources and hearing the amplified audio through headphones.

You don't mention the most important element of the budget: Time.
I would strongly recommend to the team to read about beamforming and phased arrays. These concepts play an important role in determining the specifications of the system. In turn, the specifications of the system will determine if it is viable within the given resources (time, money, expertise, other).
An important aspect of this system is timing. If you take a closer look at the phased array concept, you will notice that direction of arrival depends on time of arrival of wavefronts at each receiver (microphone). Consequently, this means that two things on the array must be fixed: Sampling Frequency and Position.
If there is jitter in the sampling frequency, then the ability of the phased array to resolve between two nearby directions is reduced. Similarly, if the microphones are not fixed (worn on body), the end result will be the same.
Therefore, one key question to ask yourselves is: "What sort of angular resolution do we want the phased array to be able to resolve?". And by the sound of it, you have to think about this in two dimensions: One, rotation around the longitudinal axis of the body and two, elevation around the lateral axis of the body. The angular resolution will tell you:

How many microphones do you need;
What Sampling frequency to run them at.

Let's park this here for a minute and approach it from the other side:
You are going for a MEMS I2S Microphone, which is not a bad idea (on its own). Adafruit's library reads the microphone with 32 bits of precision (page 11), but the manufacturer of the SPH0645 says it has 18 bits of precision and it is read at frames of 24bits (page 3). Furthermore, the I2S bus, specifies that each data line can "carry" two channels. In fact, Adafruit's website mentions that it is possible to plug two of these mics on the same I2S bus but on different "channels" for stereo operation.
So, on the hardware side you have the following things to consider:

How are you going to simultaneously trigger the sampling of $N_{mic} \ge 8$ microphones if they are on at least 4 different I2S channels?

How are you going to transfer / handle the bandwidth of $N_{mic} \times wordLength \times F_s$ bits/sec to the CPU and process it? (Where wordLength is in bits, e.g. 24bit)

Point #1 Is extremely important because of the jitter mentioned above. Adafruit's code works with polling. In other words, you set up a counter that triggers an Interrupt Service Routine and within that routine you trigger a read on the I2S bus. If you trigger a read on the first I2S bus (channels 1,2) and after you are done with them, trigger a read on the second I2S bus (channels 3,4), the values that you read from the second bus are already in the future with respect to the values of the first.
So what?!?
Let's assume that you run your I2S at 4MHz. Ignoring the bus control signals, 24 bits of samples take approximately $\frac{1}{4000000Hz} \times 24 \approx 0.000006 sec$ (time to transmit both channels). This means that if you run your system at an $F_s$ of $16kHz$ then, in one period you can read at most $\frac{\frac{1}{16000}}{0.000006} \approx 10$ channels. But this is still pushing it, because the last two channels will be out of phase to the first two channels by about 1 period of $F_s$. Is this significant? It depends on your angular resolution and the sort of phase differences that you want to measure.
Now it is the turn of point #2 to start biting. We just spend about a period of the sampling frequency in reading data in from the I2S buses attached to our microphones. Who is going to do the processing? While the ISR is executing, the CPU is not doing anything else (it is assumed of course here that once you are in the sampling ISR, you turn off all other interrupt triggers).
So, somehow we have to make up time for the processing. How much time? I don't know, it depends on what you are trying to do (beamforming!). The simplest form of a phased array sums delayed versions of the signals from the microphones. This is not costly but if you want to do it in realtime you would have to read the two angles from the user. Which means that now, you have to get the CPU to read two pots (at least) on some analog input to determine where around the user is the main lobe pointing at. So, that's more time that has to be taken OUT of the processing and the reading data in.
Hold on!, didn't we say earlier that we are supposed to store the data to "...some storage device..."? So, this means more time to write the bits to the SD card, PLUS more time to handle the FAT Filesystem (at least), to be able to write to an SD card that is then plugged to a computer and it is possible to read files off of it.
Can you see the headache? At all times you are managing timing and bandwidth. In the end, you have a DSP CPU that can do very fast operations but you have it doing chores like "read a pot", "read a mic", etc.
So, forget about MEMS Mics on I2S buses and multiplexers. If you want to do this properly, you would design a "sound board" around a codec. For example this one, 24bit 6 channel in with max $F_s$ at 192kHz and I2S connection to the CPU. Its analog Input / Output section is built to accept Line-In levels which means that you can use simple capacitive mics with simple preamps. It operates on "frames" and it has buffered mode as well, so the CPU can now process a batch of samples in one go. Of course, now, you still have to solve the problem of storage and processing multi-channel streams.
But by the time that you are into designing a sound board around a codec, you are very close to using a sound card. Check out for example the "sound card" of the gumstix, it is just this codec but already plugged into the I2S of its ARM CPU, PLUS the drivers it needs to be accessible from within the Operating System.
This is why I suggested earlier that you do away with all of this by using off the shelf items. At least, you get all the "plumbing" for free and you can focus on developing the software part.
Hope this helps.


Discontinuities in the FFT So I am taking the Fast Fourier Transform of the following function:

$$ x[n] = \displaystyle\sum\limits_{i=0}^{5} A_{i} \cos\left(\frac{\omega_{i}}{\omega_{s}} n + \phi_{i}\right) $$

Where the sampling rate $\omega_{s} = 2 \pi 8000$, the frequencies $\omega$, the phases $\phi$ and the amplitudes $A$ are:
\begin{align}
\omega &= 2 \pi \{500, 522, 610, 675, 746, 825 \}\\
\phi &= \{0,0,0,0,0,0\}\\
A &= \{1, 1, -1, -1, 1, -1\}
\end{align}

The length of the signal is $N = 2048$ samples. No window is applied before taking the FFT.

The resulting FFT spectrum should have only a real component since the imaginary component should be zero. However, when taking the FFT and plotting the spectrum I get the following:

Why is the imaginary part not zero? It seems to be a discontinuity. However, I am not sure why it exists. Why does the imaginary part of the spectrum look like this?

Your parameters aren't correct for producing a whole number of cycles for each component.

For each $i$ the value of $ \frac{\omega_{i}}{\omega_{s}} N $ has to be a multiple of $ 2 \pi $.

Hope this helps.

Ced

Followup:

Suppose your signal is

$$ x[n] = A \cos( \alpha n + \phi ) $$

If $ \phi \neq 0 $ there is a possible solution which will give you a strictly real valued DFT bins besides one with a whole number of cycles.  The condition that must be met is: $ \alpha N + 2 \phi $ must be a multiple of $ 2 \pi $.

For a whole number of cycles, your only non-zero bins will be at $ k $ and $ N-k $ where $ k $ is the number of cylces in the frame.  There is no leakage.  The bins will be real valued if $\phi = 0 $ or $\phi = \pi $.

For the alternative above, there will be leakage and it will all be real valued.

Edit:  A simpler way of stating this is that your signal has to be symmetric about the center of your interval.

Response to comment:

Sorry, I didn't see your comment until after I posted the followup.

The mathematical answer lies in equations (23) through (25) in my blog article "DFT Bin Value Formulas for Pure Real Tones".

In order to have strictly real valued bin values, the value of $ U $ must be zero.  For whole numbers of cycles, $ V $ is also zero and eq (25) becomes undefined at the relevant bins and you have to use equation (19).  The alternative solution in my followup makes $ U $ zero, but not $ V $.


Check the stability of the system, $y[n]=u[n]$ How do I check the following system
$$
y[n]=u[n]
$$
is BIBO stable or not ?

$u[n]$ is the unit step function

My Attempt:

For the BIBO stability, the necessary and sufficient condition is $$\sum_{n=-\infty}^{+\infty} |h[n]|<\infty$$.

Let $x[n]=\delta[n]$, so $y[n]=h[n]=u[n]$
$$
\sum_{n=-\infty}^{+\infty} |h[n]|<\infty=\sum_{n=0}^{+\infty} |u[n]|=1+1+1+......=\infty
$$

which proves it is an unstable system. Is it the right way to approach the problem or do I need to first write $u[n]$ in terms of the delta function ?

(Long Comment) Sorry guys, I did not get why that system created so many conceptual problems to you. For instance, consider the unforced system $y(n)=0$. If I define $\pi:\{0\}\times\mathbb Z\to\{0\}$, $\pi(y,n)=0$ then $(\{0\},\mathbb Z,\pi)$ defines a dynamical system on the state space $\{0\}$. In fact $\pi$ is continuous, and it satisfies the identity axiom ($\pi(x,0)=x$ for all $x\in\{0\}$) and the group axiom ($\pi(\pi(x,n_1),n_2) = \pi(x,n_1+n_2)$ for all $(x,n_1),(x,n_2)\in\{0\}\times\mathbb Z$).

I agree that this is a degenerate case. I guess I can extend the state space to $\mathbb R$ by defining the system in a different way. For instance, let
$$
x(n+1)=0\qquad (1)
$$
be a ordinary difference equation and define $\pi$ as the flow of (1). Then $(\mathbb R,\mathbb N,\pi)$ defines a dynamical system on $\mathbb R$. Take the forced system
\begin{align*}
x(n+1) &= 0\cdot x(n) + 0\cdot u(n)\\
y(n) &= 0\cdot x(n) + 1\cdot u(n)\qquad\qquad (2)
\end{align*}
That makes sense to me. It is a LTI system with zero relative degree. To check BIBO stability let take a reference input $\bar u:\mathbb N\to\mathbb R$ and a perturbation $\delta u:\mathbb N\to\mathbb R$. Let $u(n)=\bar u(n)+\delta u(n)$.
Let $\bar y:\mathbb N\to\mathbb R$ be the solution corresponding to the input $\bar u$ and let $y:\mathbb N\to\mathbb R$ be the one corresponding to the input $u$. Then
\begin{align*}
\|y(n)-\bar y(n)\| = \|u(n)-\bar u(n)\|=\|\delta u(n)\|
\end{align*}
So trivially, for any $\epsilon>0$ there exists $\delta_\epsilon>0$ such that
$\|\delta u(n)\|\le \delta_\epsilon\implies \|y(n)-\bar y(n)\|\le\epsilon$. You can see this also as a continuity property of the solutions. In fact, you can define the space $Y$ of all the bounded output trajectories $y:\mathbb N\to\mathbb R$ and the space $U$ of all the bounded input ones. Let endow $Y$ with the norm $\|y\|:=\sup_{n\in\mathbb N}|y(n)|$ and $U$ with the norm $\|u\|:=\sup_{n\in\mathbb N}|u(n)|$. Then the system (2) definest an operator $T:U\to Y$, and asking for BIBO stability is the same as asking for the continuity of $T$ in the sup norms.


Help me in understanding smoothing Follow up to the question here

This is a screen shot of an intermediate step in the middle of my calculation

The dotted line called MY is supposed to be a smoothed version of db pow y
db pow y is the magnitude squared represented in base 2 log scale, of the first half of a 128 point FFT
In the every next iteration, I can see that db pow y changes drastically (it is the sound input data, so it may be anything at every step) while the graph for MY changes very very slightly

my question is , since MY is a smoothed version of db pow y, shouldn't it be running close to the graph of db pow y? I know that it will definitely not be closer to db pow y as it changes drastically at each iteration and MY changes very slightly, but I still do not understand exactly how this is supposed to work

Also, at points such as A in the graph above, the dotted line seems to "ignore" sharp variations in db pow y (as I would expect, since MY is smoothed version of db pow y), but then again, at positions such as B, the dotted line is the one changing drastically

Then, could you please explain to me how this is supposed to work?

This is the smoothing equation used:

    pDashS(i) = 0.2 * pDashS(i) + 0.8 * ((db_pow_y(i)-pres_minima(i)) > diff_between_db_pow_y_and_pres_min);
    alphaDTildS(i) = 0.95 + 0.05 * pDashS(i);
    MY(i) = alphaDTildS(i) * pastMY(i) + (1- alphaDTildS(i)) * db_pow_y(i);
    pastMY(i) = MY(i);

where i ranges from 1 to 64, half of a 128 point FFT

the actual line performing the smoothing is

  MY(i) = alphaDTildS(i) * pastMY(i) + (1- alphaDTildS(i)) *
  db_pow_y(i);

which, from what I understand from answers to my previous questions, this is an example of a low pass IIR filter which performs smoothing on db_pow_y and alphaDTildS is the smoothing parameter

Please help

Thank you

(PS: code is in Matlab)

Your recursive smoothing equation will respond slowly to sudden fluctuations (that's why it's called smoothing!).

Try and change to:

alphaDTildS(i) = 0.7 + 0.3 * pDashS(i);

That should make it closer to the underlying signal - but also less "smooth"


Decomposing an Image into High Frequency and Low Frequency Parts Using Bilateral Filter How can we decompose any image (e.g. medical image, satellite image) into high- frequency and low-frequency parts by using a bilateral filter?

Since the Bilateral Filter is Spatial Invariant and Non Linear (The output is linear combination of the data, yet the weights depends on the data in Non Linear fashion) it can't be formulated as a Spatially Invariant Filter / Convolution.

So the Low Pass / High Pass concept in it classic form doesn't hold here.

Yet you can still employ the Image Pyramid concept here with your operator being the Bilateral filter and everything will work perfectly (Namely perfect reconstruction).

It is actually very effective.


How do I differentiate between signals from an unknown program? I believe this question is on topic for this site under Algorithm debugging questions, but if not, I would appreciate a suggestion about where to post this question.

I have been messing around with detecting a car key fob that was given to me from a room mate. This key fob goes to a broken car in a junk yard and he got it for free from there. I have managed to figure out how to sort out a car key signal from background noise by filtering based on the decibel level. I am very, very, new to signal processing and I want to learn more. The problem is, I cannot figure out how to make my program tell the difference between button presses and and different car keys.

I have written this program in Python 3 and I use a RTL SDR to detect the key fob. The key, while it is supposed to operate on 315 MHz, it appears to operate on 314,873 MHz. I figured this much out using Gqrx. The Key Fob is a Ford 2S4T-15K601-AA.

I hear some keys use rolling codes and I do not know if this key does. I am not interested in actually cracking any cryptography on this key unless it is necessary to detect which key buttons are pressed.

How would I differentiate between key presses and key fobs?

This is an example of the program receiving the lock button being pressed.

➜  remote git:(master) python3 remote.py
Devices: ['00000001']
Found Rafael Micro R820T tuner
[R82XX] PLL not locked!

Wanted Frequency: 315000000.0 Hz! Actual Frequency: 314873000.0 Hz!
Offset: -127000 and Squelch: -10.0

Signal!!! Decibel: -4.172264183610216
Signal: [1]
Signal!!! Decibel: -7.696191692271029
Signal: [1, 1]
Signal!!! Decibel: -3.3640835791737893
Signal: [1, 1, 1]
Signal!!! Decibel: -4.760667362463348
Signal: [1, 1, 1, 1]
Signal!!! Decibel: -7.167746085106671
Signal: [1, 1, 1, 1, 1]
Signal!!! Decibel: -3.007769020399573
Signal: [1, 1, 1, 1, 1, 1]

Just a few pointers, since a complete answer would be way too long.

The link you give for the key fob lists its FCC ID: CWTWB1U331. Using a search engine, you can find its entry in the FCC database: https://fccid.io/CWTWB1U331

The FCC database is a goldmine because it lists many of the operational parameters of the device. For example, you can verify that the frequency range you're looking at is correct.

Among other things, you can learn that the modulation is ASK, and that the line encoding is Rolling/Manchester. In particular, the "Test Report" pdf has a couple of extremely helpful figures, which will help you narrow down the transmitted frequency and the line encoding format.

For example, from Figure 2 in that document, it seems like there are 10 bits in 5 milliseconds, for a data rate of 2000 bits per second.

Another thing to note is that essentially all key fobs will transmit a different word every time the key is pressed, so the job of identifying the key will imply finding each key's patterns, which will not be easy -- otherwise, car theft would be a piece of cake.

Finally, there are several resources online for reverse-engineering wireless signals, such as remote controls and garage door openers. I recommend that you look for those.


Is there any advantage on using a low span value in a RRC filter? I know that in a RRC filter a high value of the span gives a better response, in the sense that the RRC filter response is more near to the ideal RRC filter response. However, would there be any advantage on using a small span?

Lets say, for example, that I am sending 100 BPSK symbols, a value of span = 10, would be better, in terms of the BER, than a span = 2.

Q: But would there be any advantage with span=2?

The advantage of a RRC filter with a smaller span is that it has fewer taps, so the filter requires fewer multiplies and additions per sample.  Longer filters approximate the ideal RRC response more closely but require more computation.  Shorter filters do not approximate the ideal RRC as well but require less computation.  Typically, one will choose a filter length that strikes a balance between computational demands and performance.


Image Brightness on 2D Fourier Transformation What effect does the brightness have on the frequency spectrum of 2D Fourier Transformation?

Example. Suppose we have a gray image and calculate 2D Fourier Transformation. Then we increase the brightness of the same image for 10% and repeat 2D Fourier how does the frequency spectrum change compare to the first one?

If I understand correctly only phase would change not frequency. Is that correct?

If I understand correctly only phase would change not frequency. Is that correct?

Assuming "increasing brightness by 10%" just means multiplying each pixel's value with 1.1:

No, the phase won't change either. Since the Fourier Transform is a linear operation, you just multiply every point in the spectrum by the same factor. That doesn't affect phase, since 1.1 is a real number.


Known problems with IIR DF2 when using 64-bit floating point processing? So far I was using DF1 for biquad filters for audio DSP, but it turns out that it is possible to gain significant performance gains using the new AVX/AVX512 instruction sets using DF2. Are there any relevant problems assuming we are using 64-bit floating point precision?

Depends on how close your poles are to the unit circle. For DF2, the transfer function from the input to the state variables are given by the poles, so there can be substantial gains in there. For typical audio filters, it's not unusual to see 100+ dB of gain and that would definitely be a problem for 32-bit floating point.

These enormous gains also make updating the filters dynamically very tricky since it will typically generate massive transients.

I suggest looking at transposed form 2. It's often as efficient as DF2 but has  much better behaved state variables.


Using FFT to extract a specific frequency range in Matlab? I have a signal and from that I want to extract a specific frequency range, e.g. 10Hz to 15Hz. Now I heard that I should use Fast Fourier Transformation.

I'm going to do that in Matlab, but I have no idea where to start. Is there any good example how I can do this by passing the original signal and the frequency band as parameter?

To analyze a certain frequency band, you first have to know the sampling frequency $f_\mathrm{S}$ that has been used to acquire your signal $x$. To compute the discrete Fourier transform (DFT) in Matlab:

y = fft(x);

The output of Matlab's FFT function has length $N$ and begins with frequency 0, $N$ beeing the length of $x$. The frequency range you're looking for therefore lies in the index range
$$
N\frac{f_1}{f_\mathrm{S}} + 1\ldots N\frac{f_2}{f_\mathrm{S}} +1
$$
Where $f_1$ is the lower and $f_2$ is the upper limit of your frequency range and $f_1 \leq f_2 \leq f_\mathrm{S}/2$ must hold. The addition of 1 accounts for Matlab indices beginning at 1, not at 0. Also note that in general the above expressions are non-integer numbers.


Variance of N samples of Noise What is the variance of $N$ samples of white Gaussian noise?

I have an Equation for a continuous time which includes an $\mathrm{SNR}$ term inside it.

Now when I am going to sample it, what is the value of the sampled $\mathrm{SNR}$?

$\mathrm{SNR}$ here is $\frac{A^2}{\sigma^2}$.

What is $\sigma^2$ of the sampled version??

Actually we have the following case
Variance(parameter) = 90/p^2 * SNR * T^5
where SNR = Amplitude^2/ Sigma^2.
T is the length of the signal.
Now if we want to compute the variance of the same parameter in the discrete case we can write
Variance(parameter) = 90/p^2 * SNR * (N*Ts)^5
What is the value of SNR here ??

$\sigma^2$ of the sampled version is $\sigma^2$ except for a possible scaling of units (such as volts-squared to counts-squared)

When you sample the signal the variance does not change until subsequent digital processing such as filtering is done (assuming the sampling rate and the signal itself are uncorrelated).

Consider the simple case of a sine wave; in the analog the variance is $V_p^2/2$ where $V_p$ is the peak amplitude. If we sample the sinewave for a long enough time with an incommensurate (not an integer multiple or sub-multiple of the sine wave) sampling clock, we can see without much trouble that the statistical distribution characteristics of the sine wave (histogram, variance, standard deviation etc) are unchanged.

Further, what is very interesting is if we sample a signal with a spectral power distribution that goes beyond our sample rate (and we did not implement any sort of anti-alias filter that would have reduced the overall variance), the variance of the sampled signal (again if the sampling clock is not correlated to the signal) is unchanged!  This again is very clear by imagining the histogram of the analog and digital signals; in both cases analog and digital, with enough samples from the analog signal the histograms would have the same Gaussian shape and width (standard deviation) and variance!

This is perhaps made clearer from the plot below of samples of a  Gaussian distributed waveform. Regardless of the sampling rate, as long as we sample long enough to establish a reasonable statistical distribution, we will still result in the same Gaussian shape with the same standard deviation.

It is because of this reason that I take great precautions in making spectral density estimations from digital samples, to ensure that the higher Nyquist zones have been properly filtered prior to sampling. Once digital, without knowledge of the anti-alias filtering that was done, there is no way to tell if the spectral noise density is identical to the analog at a particular frequency offset, or the power summation of multiple Nyquist zones that have aliased into the band! However (and interesting!) if I am only interested in the total power (variance) of an analog signal, regardless of frequency location of that power (spectral density), then I can sample that signal at virtually any rate (again important to know that the sampling clock and the signal are not correlated, which for AWGN as in your case they would not be) and measure the variance of the digital signal and no that it is equal to the variance of the analog signal. As might now be clear, if the digital sampling process low pass filters the analog signal (which in practice for a "white" noise signal it absolutely will!), then the digital variance will be less than the analog. To avoid this, we would make sure that the analog signal is first low pass filtered to below the analog bandwidth of the digital sampling system (which may and often is much higher than the sampling rate).


Sampling frequency of voice and audio I read somewhere that the signal of speech has a bandwidth from 300 - 3400 Hz(Why), and audio files have a bandwidth from 50 Hz to 50000 HZ(Why).

Could someone help me , why are those constant sampling frequencies ?

It's the human hearing which defines associated bandwidths. Then the sampling frequency will be chosen to be larger than twice the bandwidth for satisfying the sampling theory and still as small as possible to reduce processing costs.

Human hearing for the general sound signals is accepted to be between 20 Hz and 20 kHz. The limits depending on individual, age, gender etc.

Note that various mechanical devices, acoustical instruments or alike can produce sound waves that are well in excess of 20 kHz upper limit and include ultrasound, however, an ordinary human will not be able to hear that, that's why the upper limit is taken to be 20 kHz.

For speech used in old analog communication (telephony audio) network, the intention was to find the minimum intelligible bandwidth to reduce the technical complexity and economic cost of building that network. This bandwidth is found to be within 300-3000 Hz. However as technology improved and other means of speech communication emerged, higher bandwidth audio was also used for better quality transmission. Digital codecs, radio broadcast etc. can use 8 kHz, 16 kHz etc. bandwidths too.


Find the impulse response of an LTI system? why there is a 3/2 in the (b) problem? Does the impulse response even exist?

For the case a-)
I assume you are trained enough in DSP to see that :
$$ \mathcal{ F } \{ \delta[n-d] \} = e^{-j \omega d } $$ for all integer $d$.

Hence given a frequency response of $H(\omega) = e^{-j3\omega }$ it's apparent that the corresponding impulse response is $$h_a[n] = \delta[n-3] $$

The problem in case b-) is in the fact that it suggests a non-integer amount of shift $\alpha$ of the unit impulse $\delta[n]$ so as to obtain $\delta[n-\alpha]$ which would then have a corresponding frequency response $H_b(\omega) = e^{ -j \omega \alpha }$. But this makes no sense in the domain of discrete-time sequences which cannot be shifted by non-integer amounts which might force you
to argue that the corresponding impulse response does not even exist.

The solution requires an investigation of the relation between continuous-time and discrete-time signals through sampling as the other answer outlines. Instead here I put a shorthand result.

First observe that for any integer $d$:
$$ \delta[n-d] = \frac{ \sin(\pi(n-d)) }{\pi(n-d) } = \text{sinc}(n-d)$$ is satisfied.

The righthand side is a sampled (and therefore discrete) sinc pulse whose continuous equivalent is $$\text{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$$

Furthemore observe that the discrete-time sinc pulse $\sin(\pi n) /{ \pi n} $ is actually the impulse response of an ideal lowpass filter with a cutoff frequency of $\omega_c = \pi$ which is indeed an allpass filter with a frequency response of:
$$ H_{apf}(\omega) = 1 ~~~,~~~ \text{for all } ~~~|\omega| < \pi $$

As a consequence, the frequency response of a shifted sinc pulse
$ \frac{ \sin(\pi(n-d)) }{\pi(n-d) }$ will be:

$$H_s(\omega) = e^{-j \omega d } $$ for any $d$, integer or real, according to the time-shift property of the Fourier transform.

Hence it can be concluded by the inverse Fourier transform relation of the last statement that the impulse response $h_b[n]$ associated with the frequency response $$H_b(\omega) = e^{-j \omega (3/2) }$$ is
$$h_b[n] =  \frac{ \sin(\pi(n-3/2)) }{\pi(n-3/2) } $$

Note that $h_b[n]$ cannot be obtained from $ \frac{ \sin(\pi n) }{\pi n}$ by just shifting it $3/2$ samples which is not possible. Therefore you should obtain the samples of the impulse response $h_b[n]$ by evaluating its formula for each n.


Downsampling in DWT algorithm I was working on wavelet signal decomposition and got confused with the downsampling part of Discrete Wavelet Transform algorithm.

If we consider a signal $\{a_0[0], ... ,a_0[N-1]\}$ of length N and use formula
$$
a_1[k] = \sum_m h[m-2k]a_0[m]
$$
then for two low pass filters
$$\{h[0]=h_0,~h[1]=h_1, ~h[2]=h_2, ~h[3]=h_3\}$$
and
$$\{g[-2]=h_0,~g[-1]=h_1, ~g[0]=h_2, ~g[1]=h_3\}$$
(which are the same thing except the second filter is time-shifted version of the first one) we will get different results.

Let me clarify. After convolution of $a_0$ and $h$ we will get $N+3$ samples $\tilde a_1[k]$, $k=-1,...,N+1$ with the first non zero sample being $\tilde a_1[-1]$. So, after downsampling we will result in the following sequence
$$
\tilde a_1[0], \tilde a_1[2], \tilde a_1[4], ..., \tilde a_1[N_1]
$$
where $N_1$ is the integer part of $(N+1)/2$, with total number of samples equal to $N_1 + 1$.

Now if we convolve $g$ and $a_0$ again we will get $N+3$ samples, but this time $\hat a_1[k]$, $k=0,...,N+2$. And after downsampling we will get
$$
\hat a_1[0], \hat a_1[2], \hat a_1[4], ..., \hat a_1[N_2]
$$
where $N_2$ is the integer part of $(N+2)/2$, with total number of samples equal to $N_2 + 1$.

So as a result, if number of samples $N$ of input signal $a_0$ is even, then number of approximation coefficients for $h$ and $g$ is different (integer part of $(N+1)/2$ is less than integer part of $(N+2)/2$). More than that, second sequence consists of exactly those samples which are excluded from the first one and vice versa.

So my questions are:

1) Is this correct or do I have a mistake somewhere?

2) If this is correct, then why implementations of dwt algorithm ignore this fact? (I've tested dwt from PyWavelets package)

Wavelet analysis results in information about both frequency and time localization.  Since the two filters are slightly offset in time, if your signal varies over time, one might expect the two wavelet results to be different.


How to Emulate Sample Frequency Offset for OFDM in Matlab? I am studying the effect of sampling frequency offset on OFDM, and in order to emulate SFO, I have to slightly change the sampling rate of the signal. In order to have a +5ppm offset, I have to change the sampling rate (in Matlab) from 1 to 1.000005. I tried to use "y=resample(x,200001,200000)" in Matlab but apparently there is a limitation and I could not do it.
How can I simulate an SFO?
How can I use resample in this case?

thanks

You effectively want a interpolated waveform that is interpolated by 200,000 such that for each new sample with the 5 ppm offset you can select one additional offset to induce that time offset, for example

x[1], x[200,002], x[400,003], x[600,004]....

(Or equally one less if your sampling clock increased in frequency 5ppm).

One very simple approach to do this is to use linear interpolation to determine the value of each new sample given the adjacent two samples, or curve fitting to more adjacent samples using higher order polynomials which is what is done with detailed filter design approaches. If your waveform is sufficiently oversampled to start with, then a simple linear interpolation approach may be more than sufficient. (To assess this you could measure the error vector magnitude of the mid sample using linear interpolation compared to a resample by 2 and see if that error is less than your target SNR requirement).

With a 5ppm lower clock the new samples would be:

$$y[n]=x[n] + mod(n,200,000)(x[n+1]-x[n])/200000$$

For n as the index of the original N sample waveform 0 to N-1

For example, If the first 4 samples were

2, 5, 8, 11

First sample: 2

Second sample: $5+ 1(8-5)/200,000$

Third sample: $8+ 2(11-8)/200,000$


Small Scale Fading time varying or position varying I would like to ask how small scale fading can affect a transmission with time.
The particular system is an antenna transmitting an ASK modulated signal at 915MHz and with bandwidth between 200KHz to 500KHz.

The receiver is static (does not move) and is placed at a distance not greater than 12 meters.

The question is, if the communication consists of slots of time of milliseconds, does the fading change for each time-slot or only will change when the receiver moves?

Thank you for the help

Regards

Small scale fading results from a superposition of a large number of propagation paths between TX and RX, which can be constructive or desctructive depending on their phase offset. As the phase offset varies rapidly with the propagation distance, the phase offsets between different paths are likely to vary rapidly too, leading to a pseudo-random behaviour of the received amplitude. We refer to this as small scale fading.

Technically, this behavior requires some form of change to the propagation conditions. So again, technically, if you did an experiment in a closed chamber with noone inside and no movements, there would not be any fading but the amplitude would remain constant.

However, in reality, there is always some movement. Most of the propagation paths are due to reflections and reflecting objects can move on their own (think of wind, people, cars, etc.). So in practice, even if TX and RX are stationary, you should expect small-scale fading to be present. A way to reduce it is using directive antennas that put stronger emphasis on the line of sight component, if that is an option.


FFT time domain average vs frequency bin average I have multiple trials of physiological data. I am doing a frequency based analysis to analyze power (amplitude) in certain frequencies of interest. Is averaging multiple trials of equal length and then taking a single FFT of the averaged signal vs computing FFT for each trial and then averaging the frequency bins the same? In practice im finding this not to be the case.

Specifically, the signal naturally has a strong 1/f component and this gets emphasized if I compute the FFT of each individual trial and then average the amplitudes (real part) of each frequency bin. Are the two equivalent? is there a right way to do things? or under what principled conditions should the choice between time domain averaging vs frequency bin averaging be made?

Let me clarify.

Fourier transform does not represent the histogram of the signal.  Fourier transform is a linear transform that takes signal from time domain (complex function) to frequency domain (another complex function).  It takes a complex function to another complex function.
Fourier transform is linear as the poster above pointed out.
Phase in your samples matters as pointed out above.  If trial-by-trial data varies in phase, then you do not want to average before doing a Fourier transform, but you also do not want to average after Fourier transform.  You want to average after Fourier transform and norm.  I will elaborate below as far as exactly what needs to be done.

The main issue here is that the question is posed wrong.  It is not "should I take the Fourier transform before averaging or after averaging".  Because it makes no difference due to linearity of Fourier transform.

The correct question to ask is "should I take the amplitude of the Fourier transform before averaging or after averaging".  For this question, the answer is before.

Here are the details.

Suppose your sampled data is represented by the sequences:

$d_1=d_1[n_1],d_1[n_2],...d_1[n_N]$

$d_2=d_2[n_1],d_2[n_2],...d_2[n_N]$

$d_3=d_3[n_1],d_3[n_2],...d_3[n_N]$

...

$d_M=d_M[n_1],d_M[n_2],...d_M[n_N]$

where $d_1,...d_M$ are data from M trials and $n_1,...n_N$ are sampled time points, then:

$F_1 = \sum_{j=1}^M{|\mathcal{F}\{d_j\}|} \neq  |\mathcal{F}\{\sum_{j=1}^M{d_j}\}| = F_2$

So while the transform $\mathcal{F}$ is linear, $|\mathcal{F}|$ is not.

Furthermore, while $d_j[n_i]$ is real for all $i,j$, $\mathcal{F}\{d_j\}$ is not, but $|\mathcal{F}\{d_j\}|$ is.

As for what you should do, you should take the Fourier transform of individual trials (via FFT), get the amplitude of individual trials, and a average them together.

Finally, what is $1/f$.  $1/f$ is a short term for the frequency spectrum of "natural" signals (usually people think of images).

When people say there is a large $1/f$ component, it means that the amplitude as a function of frequency looks like $1/f$.  It's totally hand-wavy... probably coming from a biologist :p

The inverse Fourier transform of $1/f$ is some sign function, but that is useless.  It is an imaginary sign function!  Real functions generate symmetric Fourier transform.

In fact saying that the spectrum is $1/f$, tells you something about the signal, but it doesn't let you recover the signal.  All that you know is that $|\mathcal{F}\{x(t)\}| = |1/f|$.  This doesn't let you uniquely determine $x(t)$ because all the phase information is gone, and we know that the structure of a signal relies heavily on its phase.

What does $1/f$ tell you?  Simply that it contains a lot of low frequency and a little high frequency.

Just as important a question, what does averaging buy you?  and more important is how to interpret the result?  Tune in tomorrow for a more in depth discussion :p


change convolution function if signal is above threshold? Let's say I have a signal $s(t)$ and two filters $f_1(t)$ and $f_2(t)$, I also have a threshold $A$.

Now I define $a_1(t)$ as $a_1(t) = \min(s(t),A)$, then I do the convolution of $a_1* f_1$.

Now I also define $a_2(t)$ as $a_2(t) = \max(0, s(t)-A)$, (ie. whatever is above $A$) then I do the convolution of $a_2*f_2$.

The total result $a$ is $a=a_1*f_1+a_2*f_2$.

Is there a name for such a kind of filter?

Basically, I do some convolution filter on some part of the input, and some other convolution on the remaining part of the input, then add the results

This can also by generalized with more than 2 filters, or continuous thresholds and continuously parametrized convolution functions

Is there a name for such a kind of filter?

No. (Ok, there probably is, but it's not commonly used)

The point is that your filter is not a linear system, and thus has properties that make it unsuitable for most applications. Hence, I'd neither heard of someone using it, nor of someone giving it a name.


Moving average algorithm that conserves integral Suppose I have point array $y_i$ of size $N$. How to implement moving average algorithm that conserves quantity
$$
I = \sum_{i=1}^{N}y_i
$$
NOTE: I don't want time shift so I would prefer to use symmetric SMA in the middle:
$$
\hat{y}_i = \frac{1}{n}\sum_{k = i-\frac{n-1}{2}}^{i+\frac{n-1}{2}}y_k
$$
where $\hat{y}_i$  - filtered point, $n$ - window size of SMA.

It is easy to show that the above algorithm conserves $I$ "in the middle" but problem is when processing edges of a signal. So question is how to process edges of a signal in order to conserve $I$?

After several failures to construct simple method to do it I am here.

Why do you want to preserve the the sum? Is it to make the filter with the same gain everywhere?

If you want to keep the sum the same you need to pad the beginning and the end of the signal as:
$$
.... ,y_2,y_1,y_0 ---- y_0,y_1, y_2,...., y_n ---- y_n, y_{n-1},y_{n-2}...
$$
Example:
suppose that you have 6 samples

$$
y_0, y_1, y_2, y_3,y_4, y_5
$$

and you want to to do a moving average with length 5:

Special treatment at the beginning:

$$
\hat y_0   = (y_1  +  y_0 +  y_0 +  y_1 + y_2)/5,
$$
$$
\hat y_1   = (y_0  +  y_0 +  y_1 +  y_2 + y_3)/5
$$

Normal filtering:

$$
\hat y_2   = (y_0  +  y_1 +   y_2 + y_3  + y_4)/5,
$$
$$
\hat y_3   = (y_1  +  y_2 +   y_3 + y_4  + y_5)/5,
$$

Special treatment at the end:

$$
\hat y_4   = (y_2  +  y_3 +   y_4 + y_5  + y_5)/5,
$$
$$
\hat y_5   = (y_3  +  y_4 +   y_5 +  y_5 + y_4)/5,
$$


How to compare 2 MFCC (compare them if they where created from same speaker) There are multiple speakers.
Each speaker generates multiple MFCC.
If I'm getting two different MFCC, can we know if those MFCC came from 2 different speakers or same speaker ?

How can we compare them and check this ?

You can find the distance between the MFCC features for different speakers. The distance should be zero for the same speaker and larger with different speakers.
You can use different types of distance functions or even similarity functions used generally to compare between probability density functions.
Also, corresponding to each speaker you should generate a single matrix of MFCC coefficients which is composed of multiple coefficients for each frame. So, the output MFCC should be of dimension: (number of frames, number of coefficients).
Every speaker's MFCC dimension should be the same to compare them.


Orders of complimentary butterworth filters for IMU I'm currently trying to design a complimentary filter for deployment on a drone, as part of an assessed semester project. Even though processing power is not an issue (we have GHz, lots of pipelining and all the good stuff), I am still electing to use an IIR filter in the final implementation, as this is easier to represent mathematically in my opinion.

As I have read Here in this answer, certain orders of filters do not produce monotonic unity gain in output. As such, at this stage I am mostly interested in proving, in a round-about way, which orders of filters produce monotonic unity gain. Phase difference is almost inconsequential at this stage, as I am reasonably confident in not using such a high-order filter that we hit any timing constraints.

So, I have written some initial matlab code to test my application on a first order Butterworth pair, as seen below. It is worth noting that our sampling frequency is 1 Khz, meaning the nyquist frequency is 500 Hz. As such the normalized cutoff frequency of 0.25 equates to 125 Hz, although this is arbitrary at this stage. Possible problems include the summing of the 2 filters in the last line, or the need to add some integrator term to the high-pass. We also currently have some stock complimentary filters that came with the drone (all came as a big package from Quanser, their autonomous vehicle research suite thingy), but it seems lackluster to just use them and not know what they do. I will also be taking a closer look at the filters for the IMU and magnetometer in a short while after this is complete, as I think I can do a little better.

clc;
clear;
close all;
w_n=0.25; %cutoff frequency

[zL1,pL1,kL1] = butter(1,w_n,'low');
[numL1,denL1] = zp2tf(zL1,pL1,kL1);
sysL1=tf(numL1,denL1);

[zH1,pH1,kH1] = butter(1,w_n,'high');
[numH1,denH1] = zp2tf(zH1,pH1,kH1);
sysH1=tf(numH1,denH1);

sys1=sysL1+sysH1;

Making the bode plot of this transfer function gives me something completely different from what I was expecting, as also seen below.

Any help and guidance is appreciated, many thanks!

You are simply seeing the quantization of the double type. Basically you have a pole and a zero that should cancel each other out. But because of the limited precision of the double data type, it doesn't. So instead of having 0 dB you have a gain of 10^-15 dB give or take.

But in real life,  99.99999% of the time 10^-15 dB = 0 dB

For the record, I think you could use the minreal function in Matlab to cancel the almost equal pole(s) and zero(s). You would then get your unity gain and your flat bode plot with 0 dB and 0 degree phase.

I don't like simply posting links, but a reminder on how floating-point works should help you.

https://en.wikipedia.org/wiki/Floating-point_arithmetic


rectangular function in matlab I am trying to program this function: $ g(t)= \operatorname{rect}(t/T) \,f(t)$

I take: $t=0:1:20$ and I know how to program $f(t)$.

But I have a problem with the rectangular function. How to program this function in matlab?

function y = rect(x);

    x = abs(x);

    if x == 0.5
        y = 0.5;
    else
        if  x < 0.5
            y = 1;
        else
            y = 0;
        end
    end

end

put this in a file and name it "rect.m".  make sure it's on your MATLAB search path.


Relationship between order of filter and amplitude/phase response I am learning about digital filters, for example a lowpass butterworth filter.

In many languages, the implmentation requires specification of the order of the filter, which I believe is the number of terms in the closed form equation that generates each filtered sample, given the original samples and previously generated samples.

What I haven't been able to find (at a sufficiently ELI5 level) is how the order affects the frequency and phase response of the filter.

Does a higher order simply mean a better approximation to the optimal filter with (cutoff, stopband) ? Do higher orders mean different or less severe phase differences?

Finally, I'd assume that it affects performance somewhat as well.

Does a higher order simply mean a better approximation to the optimal filter with (cutoff, stopband) ?

That depends on how you define "optimal". In general a high order will give you a steeper rolloff

  Do higher orders mean different or less severe phase differences?

Higher order will result in more phase distortion

  Finally, I'd assume that it affects performance somewhat as well.

That depends on how you define "performance". The "best" filter always depends on the application.

In general higher IIR filter order will give you more steepness in the frequency but also more time smear and phase distortion. Any beginner class on digital filter should cover this. For example, take a look at https://ccrma.stanford.edu/~jos/filters/


once again,confusion between phase and group delay Despite,reading up on this I'm unable to find crystal clear clarity on the differences b/w phase and grpdelay functions and interpretation of grpdelay's output. I have designed a low pass butterworth filter having a cut off frequency of 0.04Hz. There is a delay in my filtered output compared to my original signal. I wish to know this delay. I used the grpdelay function in matlab (grpdelay(b,a,N/2) ). And this is the output i got..
 I read that grpdelay is supposed to give the phase responses as a function of frequencies. In my image, the cursor points to the cut-off freq. How exactly do I know the delay in my filtered output through this? Is it the corresponding value on the Y-axis?

Plotting the phasedelay(b,a), gave this as the output..Here,the phasedelay value at cut-off is different.

My question is, to figure out how much the filter has delayed the signal,which function do I use? And how do I interpret the output exactly, of that function?

Here is my original signal in blue and the filtered signal which has the delay in red .

If you are looking for a frequency-independent delay applied to any given input signal by the filter (apart from amplifying and attenuating certain frequency components), then you won't be able to find it because there is no such delay. As you can see in your plots, group delay and phase delay are generally frequency dependent. Furthermore, for general input signals, the terms group delay and phase delay have no meaning.

There are two cases for which the terms group delay and phase delay have a clear meaning. The first is for narrow-band input signals

$$x(t)=a(t) \cos(\omega_0 t)$$

where $a(t)$ is a  low-pass signal. If the system has an approximately constant amplitude response in the frequency range of $x(t)$ (i.e., around $\omega_0$) which is $\left|H(j \omega_0)\right|$, and if its phase is approximately linear in this frequency range, then it can be shown that the output signal is approximately given by

$$y(t) \approx \left| H(j \omega_0) \right| a(t-\tau_g(\omega_0)) \cos\left( \omega_0(t-\tau_\phi(\omega_0)) \right)\tag{1}$$

where $\tau_g(\omega_0)$ is the system's group delay evaluated at the carrier frequency $\omega_0$, and $\tau_\phi(\omega_0)$ is the phase delay at $\omega_0$. So in this case the group delay is the delay of the envelope, whereas the phase delay equals the delay of the carrier. A special example of such a narrow band input signal is a pure sinusoid. From (1) it is clear that the delay of a pure sinusoid when passing through a linear time-invariant (LTI) filter equals the phase delay evaluated at the sinusoid's frequency (and not the group delay). Note that in the case of a pure sinusoid (i.e. $a(t)=\mbox{const}$), Equation (1) becomes exact for any LTI system.

The other case where these terms are meaningful is when the system has exactly linear phase. This is only possible for FIR systems. For FIR filters with even symmetry of the impulse response, the group delay and the phase delay are equal and constant. In this case, they simply specify the delay of the input signal. This is possible because both group delay and phase delay are constant, i.e. frequency independent.

From your question, I believe that you are looking for such a delay. However, as explained above, such delay only exists for linear-phase FIR systems, not for the Butterworth filter you are considering.


A 3x3 image kernel takes 9 multiplications and 8 or 9 additions for each pixel? I doubt because you have 9 multiplications and you are adding 9 numbers together, when using a 3x3 convolution matrix.

Adding 9 numbers: 1+2+3+4+5+6+7+8+9

We are performing 8 additions  as we repeat the operation of 'adding' 8 times (8 times a '+' sign).

We can say there are 9 additions as well, if you would accept a definition like: 0+1+2+3+4+5+6+7+8+9

But what does make more sense for computers: 8 adds or 9 adds?

As you have 9 terms to add together, your computer will perform 8 additions as each addition needs two inputs. So you have 9 multiplications and 8 additions.


What does the normalization step of the Haar wavelet transform represent? When you perform the Haar wavelet transform, you take the sums and differences, then at each stage, you multiply the entire signal by $\small\sqrt2$.

When taking the inverse transform, you multiply the signal by $\frac{1}{\sqrt2}$ for each iteration.

What does this "normalization" really represent?

As I understand it, the normalization is because the Haar wavelet conserves energy of the signal. In that, when you take signal from one domain to another, you aren't supposed to add energy to it, (although conceivably you might lose energy).

The normalization is just a way to ensure that the energy of your Haar-transformed signal in the Haar-domain has the exact same energy as your signal in the original domain.

Intuitively speaking, Haar, Fourier, etc, are all just basis transformations, which intuitively just mean that you are looking at the signal in a different way, (technically, through a different set of bases). Therefore, if all you are doing is looking at a signal differently, its energy cannot/should not change.


Applying Superposition Property : $x^{2}[n]$ vs $x[n^{2}]$ If we consider the mapping $\mathcal{H} : x[n]\mapsto y[n]$ and define the following output signal $y_{1}[n]:=\mathcal{H}\{x[n]\}:=x^{2}[n]$, then one can easily verify that such system is non-linear for it rejects the superposition property. But I am having a hard time applying the superposition property to $y_{2}[n]:=\mathcal{H}\{x[n]\}:=x[n^{2}]$ because its confuses with me with $y_{1}[n]$. I hope someone can offer any assitance I would be very much grateful.

Given the system I/O definition:
$$y[n] = \mathcal{H}\{x[n]\} = x[n^2] \tag{1} $$
you can easily show that it's a linear (but time-varying) system.
Following the standard procedure
let
$$y_1[n] = \mathcal{H}\{x_1[n]\} = x_1[n^2] \tag{2.1}$$ and
$$y_2[n] = \mathcal{H}\{x_2[n]\} = x_2[n^2] \tag{2.2}$$
then define
$$x_3[n] = a x_1[n] + b x_2[n] \tag{3}$$
and
$$
\begin{align}
y_3[n] &= \mathcal{H}\{x_3[n]\} \tag{4}\\\\
&= x_3[n^2] \tag{5}\\\\
&= a x_1[n^2] + b x_2[n^2]\tag{6}\\\\
&= a y_1[n] + b y_2[n]\tag{7}\\
\end{align}
$$
where Eqs(4)-(5) follow Eq.(1), Eq.(6) follows Eq.(3), and Eq.(7) follows Eq.(2).
Equation (7) indicates that the system is linear.


Demodulating AM signal from carrier signal on Matlab I'm trying to demodulate a voice signal (bandwidth approximately 5 kHz) from a carrier signal (100kHz).

  [y,Fs]=audioread('Signal1.wav');

  h=10*sinc(10*(-10:10));

  y1=filter(h,1,y);

But when I compare the signal with the "ground truth" signal provided, the size of the signal is 1/20th of the one I made. Did I miss some steps in the process? Thanks.

The reason of different signal lengths is due to their different sampling rates.
The original baseband analog audio signal $m(t)$ has a $5$ kHz bandwidth, and should be sampled at $F_{bb} \geq 10$ kHz, to represent it (without aliasing) in the discrete-time form $m[n] = m(n / F_{bb})$ in Matlab.
Assuming you want to use the lowpass sampling theorem, then the amplitude modulated analog signal $x(t) = A_c m(t) \cos(2\pi f_c t + \phi)$ has a carrier frequency of $f_c = 100$ kHz (its bandwidth is $100+5$ kHz) and has to be sampled by at least $F_{am} = \frac{1}{T_{am}} = 210$ kHz to represent it as $x[n] = A_c m_{am}[n] \cos(2 \pi f_c n T_{am} + \phi)$.
Note that in the latter case, the message signal $m(t)$ is oversampled into $m_{am}[n]$ by about $F_{am}/F_{bb} = 210/10 =  21$ for getting the sample by sample product. Hence the number of samples in the AM modulated (and demodulated) signal $m_{am}[n]$ which represents the audio information will be about $21$ times more than that of the original baseband audio representation in $m[n]$.
Typically you would downsample the demodulated audio sequence in order to reduce its sampling rate (and sample count).


Block LMS with overlapping blocks In the Block LMS algorithm, the input is partitioned into nonoverlapping blocks of size $L$ and the filter coefficients are updated once every $L$ samples. Would convergence improve if we used more frequent filter coefficient updates by using overlapping blocks? Having more "descent steps" per unit time would seem to imply faster convergence, but the lack of mention of this in any literature suggests otherwise.

There is no hard rule regarding convergence speed of the block-LMS vs sample-by-sample LMS. It really depends on the scenario.
On top of my head is the following two (stationary) scenarios:

A very noisy scenario, where a single estimate of the gradient is not enough. In this case, the block-LMS has better gradient estimates and would usually result in faster convergence.
A clean scenario. In this one, the sample-by-sample LMS will yield a gradient estimate that is almost as good as the block-based LMS. Then, the sample-by-sample LMS will show a better convergence behavior.

Other examples can be constructed using stationary/non-stationary distinction.


Clarification on Discrete-time convolution and signal shifting What are $y_1[n]=x_1[n]*h[n]$ and $y_2[n]=x_2[n]*h[n]$?
  $$x_1[n]=(0.1)^nu[n],\quad x_2[n]=(0.2)^n,\quad h[n]=(0.3)^nu[n]$$

I am confused about how to calculate $y_{1}[n]$. The formula for discrete convolution is:
 $$\sum_{m=-\infty}^{\infty}x[m]h[n-m] = \sum_{m=-\infty}^{\infty}h[m]x[n-m]$$.

My attempt:
I chose to shift $x_{1}[n]$ along m-axis:
$$y_{1}[n] = \sum_{m=0}^{n} (0.3)^{m}(0.1)^{n-m}$$
$$y_{1}[n] =(0.1)^{n}\sum_{m=0}^{n}(\frac{0.3}{0.1})^{m}$$
$$y_{1}[n] = (0.1)^{n}\sum_{m=0}^{n}(3)^{m}$$

This gives me $y_{1}[n] = \infty$. However, when I shift $h[n]$ along the m-axis I get $ y_{1}[n]$ as a finite value. My understanding of discrete time convolution was that shifting either $x[n]$ or $h[n]$ would result in the same answer.

You're right that both results should be identical. The mistake lies in your conclusion that $y_1[n]=\infty$ in the first case, that's not true. You get

$$y_1[n]=(0.1)^nu[n]\sum_{m=0}^n3^m=(0.1)^nu[n]\frac{3^{n+1}-1}{3-1}=\ldots\tag{1}$$

Note that you must also include the unit step $u[n]$ in $(1)$ because the response is zero for $n<0$.


How to plot magnitude response using $\tt freqz$ in MATLAB if I have a transfer function without a denominator So I have the transfer function:

$$H[z] = 1 + \sqrt{2}z^{-1} + z^{-2} $$

And I have to evaluate $H(e^{j\omega})$ for $\omega= 0, \pi/4, \pi/4 \ldots$

I have done the calculations manually using Euler's formula, but now the assignment is asking me to compare these plots with the plots using freqz in MATLAB. I can't seem to find instructions on how I can do that with this type of transfer function.

You simply specify a = 1 (because the denominator equals $1$). So you get

b = [1,sqrt(2),1];
a = 1;
N = 512;
[H,w] = freqz(b,a,N);

You can compare this to the analytical solution:

H2 = 1 + sqrt(2)*exp(-1i*w) + exp(-1i*2*w);
max(abs(H2-H))    % 8.0825e-16


dB/hz to dB Conversion I'm have a PSD graph of dB/hz vs. frequency in hz. I want to convert the dB/hz unit to dB, but I can't find any way of doing it. Help would be much appreciated!

Recall from Physics class that density is, mass per unit volume. As an example, the approximate density of water is 1 g/cm3, so 100 cubic centimetres of water weighs 100 grams.

Power Spectral Density (PSD) is simply average noise power per unit of bandwidth. Power is expressed in units of watts. In telephony measurements, power is commonly expressed in logarithmic terms with respect to a milliwatt, i.e., dBm. The unit (volume) of bandwidth is
1 Hz. So PSD may be expressed as dBm/Hz, pronounced, “dBm per Hz.”

PSD is to Total Power as Density is to Weight. All transmission systems have finite power and bandwidth. Let’s consider an ideal system with a flat PSD signature.

At each frequency between f1 and f2 Hz, Let PSD to be a flat -10 dBm/Hz, or 0.1 mW/Hz.

Then

PSD (dB) = PSD(dbm/Hz) * Bandwidth(Hz)

where bandwidth can be defined as f2 –f1.

If you want to calculate the total power there are more info available in the reference.

reference: PSD techniques by Peter Walsh, NCE


Identify the system from input and output I am using Excel. If I have an input and an output, and I know that a recursive band-pass filter produced the output, how can I can find the coefficients of that filter?

You can use an equation error method and get an estimation of the coefficients by solving an overdetermined system of linear equations. The transfer function of an IIR filter is given by

$$H(z)=\frac{B(z)}{A(z)}=\frac{\sum_{m=0}^{M}b[m]z^{-m}}{1+\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

From (1), for given input and output sequences $x[n]$ and $y[n]$, respectively, the following equation must be satisfied:

$$(a*y)[n]=(x*b)[n]\tag{2}$$

where $*$ denotes convolution. Equation (2) can be written in matrix/vector form:

$$\mathbf{y}=\mathbf{Xb}-\mathbf{Ya}\tag{3}$$

where $\mathbf{y}$ is the vector of the output signal, $\mathbf{b}$ and $\mathbf{a}$ are the unknown coefficient vectors (excluding $a[0]=1$), and $\mathbf{X}$ and $\mathbf{Y}$ are the convolution matrices consisting of the values of the input and output sequences, respectively. Equation (3) can be rewritten as

$$\mathbf{y}=\mathbf{[-Y,X]\cdot[a,b]^T}\tag{4}$$

where $^T$ denotes transposition. The overdetermined system (4) can be solved in a least squares sense using standard methods, such as the one implemented in Matlab or Octave.

Note that in any case you need to estimate/choose the orders $M$ and $N$ of the numerator and denominator polynomials in (1). For most standard filters you have $M=N$.

This little Matlab/Octave script shows how it works:

x=randn(100,1);       % some input signal
a=[1 .3 -.4 .2 .1];   % denominator coefficients
b=[2 3 4 3 2];        % numerator coefficients
y=filter(b,a,x);      % filter output
M=3;                  % guessed numerator order
N=3;                  % guessed denominator order

% build matrices and solve overdetermined linear system
A = flipud(hankel(flipud(y(:)))); A1 = A(:,1); A = A(:,2:N+1);
B = flipud(hankel(flipud(x(:)))); B = B(:,1:M+1);
C = [-A,B];
c = C\A1;
a1 = [1;c(1:N)];
b1 = c(N+1:M+N+1);

% compare outputs of original and estimated filter
y1=filter(b1,a1,x);
n=1:100;
plot(n,y,n,y1)


How to make signal smooth using filter? I have measured the output of a mechanical system and obtained the signal shown at the attached image. I want to make it smoother by removing the chattering (I think it is called chattering). In order to achieve this I used a filter with transfer function:

$\ T(s) = \frac{30}{s+30} $

I used this filter to filter the measured data and the result was the the signal at the attached image. I would like to further process the signal to make it as better as possible. Can I achieve what I want by using another type of filter or is there any other process to follow ?

From what I can glance :

The chattering seems to be high-frequency compared to your signal of interest, it should not be hard to filter this chattering noise. You simply need to identify the frequency band of your signal and the frequency band of this noise.

Could your perform an FFT to analyze the frequencies of your noise? Then design a filter that will not attenuate your signal of interest while attenuating the noise to an acceptable level.

Is it real-time processing or post-processing?


How does this FLL work? I'm trying to get my head around how band-edge frequency locked loops work. In particular, a diagram from this paper:

Specifically, what does the Tanh() function and the conjugate product of the two filters accomplish?

Explanations and references would be much appreciated.

The conjugate product of the two filters form a frequency discriminator with an error output that is proportional to a frequency offset at the input.

This is due to the relationship that the cross correlation in one domain is the conjugate product in the other domain (you can skip past the derivation unless interested):

$$c(\tau) = XCORR[x(t)y(t)]$$
$$c(\tau) = \int_{t=-\infty}^{\infty}x^*(t)y(t+\tau)dt$$
$$C(\omega) = \mathscr F\{y(t)\}= \int_{\tau=-\infty}^{\infty}\bigg[\int_{t=-\infty}^{\infty}x^*(t)y(t+\tau)dt\bigg]e^{-j\omega \tau}d\tau$$
$$=\int_{t=-\infty}^{\infty}x^*(t)\bigg[\int_{\tau=-\infty}^{\infty}y(t+\tau)e^{-j\omega \tau}d\tau\bigg]dt$$

From the shift property of the Fourier Transform, the terms inside the brackets are:

$$\int_{\tau=-\infty}^{\infty}y(t+\tau)e^{-j\omega \tau}d\tau = Y(\omega)e^{j\omega t}$$

Therefore:
$$C(\omega) =\int_{t=-\infty}^{\infty}x^*(t)\bigg[Y(\omega)e^{j\omega t}\bigg]dt$$

$$ = Y(\omega)\int_{t=-\infty}^{\infty}x^*(t)e^{j\omega t}dt  = Y(\omega)\bigg[\int_{t=-\infty}^{\infty}x(t)e^{-j\omega t}dt\bigg]^*$$

$$ =Y(\omega)X^*(\omega) $$

The above shows that the cross correlation of two functions x(t) and y(t) in the time domain is equivalent to the conjugate multiplication of their Fourier Transforms (in the frequency domain). Similarly, and this is the important point to explaining the operation of the frequency discriminator, is that the converse is also true: the cross correlation of two functions in the frequency domain is equivalent to the conjugate multiplication of the transform of the same two functions in the time domain.

Consider the following as example spectrums of the "Matched Filter" and "Frequency Derivative Matched Filter" in the block diagram given by the OP. Notably the Frequency Derivative Matched Filter is a Band-Edge Filter.

When there is no frequency error, the input spectrum would be aligned with both of these resulting in the following spectrums at the outputs of these filters. Notably, observe that convolving these two output spectrums will result in a zero result specifically due to the positive and negative balance in the Frequency Derivative Matched Filter Spectrum!

When there is a positive frequency error in the input signal as depicted in the graphic below, the input spectrum would be offset by the frequency error resulting in the following spectrums at the output of the two filters. In this case specifically, convolving the two output spectrums would result in a negative value corresponding to the frequency error.  Similarly a negative frequency error would result in a positive output value after the convolution in frequency of the two spectrums.

The above shows how the convolution of a Matched Filter Spectrum and a Frequency Derivative Matched Filter Spectrum results in an output value that is the error term for a frequency discriminator. Additionally we showed that convolution in frequency is equivalent to the complex conjugate multiplication in time. Thus the complex conjugate multiplication of the two filters in the block diagram form a frequency discriminator as the error signal into the loop.

The frequency discrimintor would operate as described without the scaling by the Signal to Noise Ratio (SNR) and Tanh() operation.

The additional scaling by the SNR, and the Tanh() operation are done to perform maximum likelihood frequency estimation (MLE). This shows the ideal implementation (in a maximum likelihood sense) but to my understanding the SNR scaling is rarely done in practice (but could be if an SNR estimation is performed as part of the receiver operation), and the Tanh() operation is typically approximated with a piecewise linear operation to be the signal for the small signal case when the signal is less than ±1 or limited to the sign of the signal (+1/-1) otherwise.  So basically this path in the block diagram is hard limited for better performance with a limit that is dependent on the SNR of the received signal (a SNR dependent hard limiter is a typical operation within MLE implementations).

I believe that fred harris (fh) points out an advantage to this approach is how the derivative matched filter response is concentrated at the band edge where as we see through the convolution operation most of the information regarding frequency offset resides (when using a complex conjugate multiplier as a discriminator). The filter naturally blanks out the center portion of the spectrum which is not contributing further the our frequency error answer but would be contributing noise otherwise. Other advantages to this approach are that the symbol timing clock (timing recovery) can also be derived. I have not studied it personally in detail nor implemented this approach, but I believe fh makes the case that this implementation results in much lower pattern noise versus other carrier recovery (and timing recovery) approaches.


camera motion detection I am working on camera motion detection.How do I decouple translation and rotation? I have obtained rotation by using sensors. I understand we can do something as simple as X'=RX+T , where R-rotation, T-translation and X,X' are the points involved. These points are 3D points(either world or camera coordinate system). How do i determine these 3D points by analyzing an image?

Is there any other way to obtain translation? (some Image Processing technique may be?)

Any links for literature would be helpful.

If X are known 3D points, then the translation is given simply by:

$$T=\frac{\sum_{i=1}^{n}{X'_{i}-RX_{i}}}{n}$$

In case of points being in projected coordinates, then more work have to be done. You can express the 3D points in homogenous coordinates:

$$X=(x, y, z, w)$$

where $w$ is 1 by default.

Then both the translation and rotation can be joined via matrix multiplication:

$$X'=TRX=\begin{bmatrix}1 & 0 & 0 & t_{x} \\ 0 & 1 & 0 & t_{y} \\ 0 & 0 & 1 & t_{z} \\ 0 & 0 & 0 & 1 \end{bmatrix}RX$$

Furthermore, you would be able to add camrea projection matrix in the equation:

$$X'=TRPX$$

Now each point correspondence ($X\leftrightarrow X')$ will give you 4 equations. The minimum number of point correspondences depends on the number of unknowns.

There are usually more than enough points and the problem is overdetermined. The above equation can then be solved in a least-squares fashion.

Another approach is to first perform camera calibration to determine matrix $P$. Then the problem will be simplified greatly.

The best reference I know is (what else) H&Z book. Chapter 7 is all about computing the camera matrix.

Estimation methods are described as well.


Convert old computer cassette signal into binary Completely new to this but keen to learn :)
I have an old BBC Micro computer cassette which I've recorded into unsigned 8-bit pcm wav file format. The file comprises tracks separated by virtual silence and each track is split into a series of blocks separated by a carrier tone at 2400hz. The blocks are a series of 1200hz and 2400hz signals representing binary data. I think this is frequency modulation? (as you can tell, I'm not very educated in these things)

Where would I start if I wanted to convert this into binary?

First of all, you'd go on an information scavenger hunt through the internet on that cassette's data format: http://beebwiki.mdfs.net/Acorn_cassette_format tells us you're quite spot on:

Signal
The signal may be in one of three states; zero, one or no carrier.
Breaks in the carrier are detected and used to reset the cassette
loading firmware.
Zero is represented by a sinusoidal wave at 1200 Hz. The exact
frequency is nominally 16 000 000 / 13 312, or 1201.9 Hz, but tape
decks vary in speed so small differences are tolerated. One cycle (at
1200 baud) or four cycles (at 300 baud) represent one zero bit.
Binary one is represented by a sinusoidal wave at 2400 Hz. The nominal
frequency is closer to 2403.8 Hz. Two cycles (at 1200 baud) or eight
cycles (at 300 baud) represent a one bit. An odd number of 2400 Hz
cycles can and does occur.
To allow the tape deck circuitry to settle, each data stream is
preceded by 5.1 seconds of 2400 Hz leader tone. This is reduced to 1.1
seconds if the recording computer has paused in the middle of a file,
or 0.9 seconds between data blocks recorded in one go.
At the end of the stream is a 5.3 second, 2400 Hz trailer tone. This
is reduced to 0.2 seconds when pausing in the middle of a file (giving
at least 1.3 seconds' delay between data blocks.) The timings are
derived from VSYNC interrupts so they vary between recordings.

So, yes, this is Frequency Shift Keying (FSK) and can be demodulated as such.
Considering modern PCs'/laptops'/phones'/fridges'/toaster ovens' computational power is basically huge compared to the data that fits on a cassette, we don't need to be very considerate on how efficiently we demodulate this, at all.
So, you'd go and

design two bandpass filters: one that's centered around 1201.9 Hz, the other one at 2403.8 Hz,
take their outputs, convert them to their absolute value, then
weight one of them so that they have the same amplitude when excited,
subtract them and finally
make a 0/1 decision based on the sign of the result to get the bit.

Sounds a bit taunting, but looks roughly like this in GNU Radio 3.9:

(You can ignore everything that says "QT" in its name, it's just for interaction / vizualization. Once you've figured out the right amplitude so that 0 and 1 bits come out the same (and that's really just necessary to minimize the probability of bit errors), you'd just multiply with that fixed value instead of having the Qt range.)
You'll have to add some discarding for the start and stop bits surrounding each byte, and you'll need to keep synchronized to these – I think the most elegant (and easy) way would be to write a quick few lines of python code in an "embedded Python Block" in GNU Radio Companion to detect the start of the start bit, then pass on through the 8 data bits, ignore the stop bit and go waiting for the next start bit. Pretty doable!


What is the impulse response for 2D functions I study image processing and I am new to it. so my question is for example what is the impulse response for:

$$ y(m,n) = 3\times x(m,n)+9 $$
$$ y(m,n) = m^2n^2 \times x(m,n) $$
where the $y(m,n)$ is output and $x(m,n)$ is input.
thanks.

As pointed out by Jason R. in his comment, you can easily apply an impulse $\delta(m,n)$ (with $\delta(m,n)=1$ for $m=n=0$, and zero otherwise) to the systems and see what the output is. However, for those two examples the response to an impulse does not say much about the general behavior of the system.

There is an important class of systems called Linear Shift Invariant (LSI) systems which are completely characterized by their impulse responses. Linearity means that if $y_1(m,n)$ is the response to an input $x_1(m,n)$ and $y_2(m,n)$ is the response to input $x_2(m,n)$, then the response to the input $ax_1(m,n)+bx_2(m,n)$ is given by $ay_1(m,n)+by_2(m,n)$. If you check this condition you will find out that your first system is not linear.

Shift invariance means that if $y(m,n)$ is the response to $x(m,n)$ then $y(m-k,n-l)$ is the response to the input $x(m-k,n-l)$. Your second system does not satisfy this condition.

So none of the two systems is LSI and therefore neither of them can be characterized by their impulse response.


How to fit an ellipse to 2D data I'd like to find the "best" fit of an ellipse to contiguous, possibly concave shapes such as:

What have I tried?

I thought that one could assign the direction of the major and minor axes $\vec a, \vec b$ of the ellipse by mapping the pixel values to coordinates, mean subtracting, and saving the largest two eigenvectors from a PCA. The seems to work rather well in finding the direction:

My problem is determining the length of these two vectors. For now, I've used $\sqrt \lambda_1$ $\sqrt \lambda_2$, from the eigenvalues of the PCA. This seems to underestimate the length. How can I determine $|\vec a|$, and $|\vec b|$ or alternatively, best fit an ellipse to these shapes?

Following up on my deleted answer... If you take a filled ellipse and project all the points onto the $x$ axis, more points will be projected near the origin than on the extremes, in a circular-shaped distribution. Not a gaussian distribution, and not the uniform distribution I mentioned in the 1-D analogy in my deleted answer. The resulting distribution actually has pdf $p(x) = \sqrt{(1 - (\frac{x}{r})^2)}$, and from there you can compute that the standard deviation is $\frac{r}{2}$.

Thus, if data is uniformly distributed in the interior of an ellipse of radii $a, b$ (whose axes are the $x$ and $y$ axes), the standard deviation of the $x$ coordinate is $\frac{a}{2}$ and of the $y$ coordinate is $\frac{b}{2}$. So the correction factor you need to use is simply 2.

Here is a worked example in python recovering the center (translation matrix), rotation matrix, and radii of an ellipse from points randomly sampled from its interior:

import numpy

# Generate some points distributed uniformely along an ellipse
x = 2 * (numpy.random.rand(10000, 1) - 0.5)
y = 2 * (numpy.random.rand(10000, 1) - 0.5)
d = (x / 0.5) ** 2 + (y / 0.25) ** 2
inside = numpy.where(d < 1)[0]
x = x[inside]
y = y[inside]
data = numpy.hstack((x, y)).T

# Now rotate by 0.5 rad and translate it to (4, -8)
angle = 0.5
rotation = numpy.array([
    [numpy.cos(0.4), -numpy.sin(0.4)],
    [numpy.sin(0.4), numpy.cos(0.4)]])

data = numpy.dot(rotation, data)
data[0, :] += 4
data[1, :] -= 8

# Step 1: center the data to get the translation vector.
print 'Translation', data.mean(axis=1)
data -= numpy.reshape(data.mean(axis=1), (2, 1))

# Step 2: perform PCA to find rotation matrix.
scatter = numpy.dot(data, data.T)
eigenvalues, transform = numpy.linalg.eig(scatter)
print 'Rotation matrix', transform

# Step 3: Rotate back the data and compute radii.
# You can also get the radii from smaller to bigger
# with 2 / numpy.sqrt(eigenvalues)
rotated = numpy.dot(numpy.linalg.inv(transform), data)
print 'Radii', 2 * rotated.std(axis=1)