Block diagram of M-PSK modulation and demodulation

So this seems to be a general block diagram for these types of modulations.

No, this is a block diagram of a IQ upconverter with some unspecified digital data modulator on the input.
You need to mentally divide the modulation from the tools you use to implement it. There's other architectures that can produce the same signal (and they're not uncommon).

But the block diagram shows that the modulator uses low-pass filters (LPF).

Yeah, that block diagram is admittedly not in a shape that I find good: it omits where the transition from digital to analog happens, where the transition from bits to symbols happens, and a few other things.
But yeah, the author of that block diagram probably meant these LPFs to signify pulse shaping filters. We can't be sure from the block diagram alone.

From what I've read, these modulations mostly use raised cosine/root raised cosine filters.

Again, no, please make sure you're not confusing the constellation (BPSK, QPSK, …) with the system, or with the pulse shaping.
RC and RRC (and basically, all pulse-shaping filters) are low-passes, so "LPF" covers these quite well, and RC and RRC are far from the only choice there, even if they are a "textbook favorite". They are popular in practice, too, indeed, but you got to realize that your textbook teaches you about the untruncated RRC and its matched filter (which happens to be itself) and their combination RC, but in reality, you have truncated versions combined with channel impulse responses and equalizers, so, make sure you don't learn that there's a "RRC in the block diagram at that position", but to understand that it's a pulse-shaping filter, and RRC is a common choice, for mathematical reasons.
A lot of your argument is "usually it's done like this and that", and reality is that there's reasons, and just knowing what is often done isn't taking you far!

Why do most of the block diagrams show blocks with this particular filter? Is it simply because it's a generic diagram?

Again, the block diagram you've shown is not that of a PSK transmitter, but that of a input data-to-bandpass signal direct converter, and not a great one at that. Of course it doesn't show your favorite pulse-shaping filter, as that choice isn't applicable to all systems covered by this, so you're absolutely right.

I have this block diagram:

That is the block diagram of a receiver architecture that is honestly doing a different kind of abstraction than your transmitter diagram, and honestly, again omits what these LPFs are for (here, they have a very important different role, namely to remove the upconversion result and obtain the baseband signal!), omits that there's some digitization going on... Again, I don't think this is a great diagram.
Also, this diagram, while underspecifying a lot of things, makes a distinct choice in terms of synchronization: The carrier recovery is done purely, it seems, on the RF signal, and without any baseband feedback, and then symbol timing is estimated on that, there's no phase recovery...
honestly, I think you need a better textbook that explains what does which in a simple circuit, instead of giving you a diagram that shows something particular and leaving it up to you, with incomplete information, to match up your knowledge on modulations with that.

Will it work for any M-PSK modulation (BPSK, QPSK, 8-PSK, 16-PSK etc.)?

No, this thing has no chance to correct phase due to the purely feed-forward nature of the carrier recovery (can't do a PLL without a loop). As such, it makes little sense to even do this in IQ.
So, doesn't work for complex constellations and certainly not for PSKs.

Or is it only valid for QPSK?

Not valid for QPSK.

Won't the block diagram look different for e.g. BPSK modulation?

Any sensible block diagram for any PSK looks different.

Will the general block diagram for M-PSK modulation look the same? And if not, what would such a diagram look like?

Your real question is "how does a PSK receiver look like" and the answer is: there's multiple architectures you can apply.
You're clearly learning complex baseband / direct conversion techniques. So, please, use your textbook (or a better one) to learn that. A typical downconverting mixer in a modern system looks something like this:

(the low-pass filters in here have a single purpose: they filter away the mixing products of the oscillator frequency $f_{LO}$ and the carrier frequency $f_c$ at $f_{LO}+f_c$ and higher intermodulation products, and limit the bandwidth of the signal so that the ADC can sample it alias-free.)
And all the frequency, phase, timing synchronization, the matched filtering, and the decision is done on the digital side, and will look differently for different systems, and depend on more factors than whether you choose QPSK or 8-PSK or 1024-QAM or 2-FSK or GMFSK or whatever.
I think the key point I want to make here is: very few use cases¹ build an analog signal chain that is specific to the modulation used in the signal in the last 20 years. In that, the diagrams you've shown are simply from a time "long past", and give you irrelevant or even obsolete detail in some aspects, and omit relevant things in others.
Learn about what complex baseband is, and what these mixers do, and the rest will become more easily understandable to you. Watch out that you don't conflate the terms of modulation, modulator, and the individual implementation of it. Generally, a good textbook's didactic approach (reading it linearly from front to back) is much faster at helping you understand what you need to do than trying to match incomplete understanding of a few things to multiple resources.
I think you're more of a practically-minded nature and it's probably a stretch assuming you'd want to learn lots of theory rather than learning the common tools of the trade and why and how they work. So, Robin Getz' and Travis Collins' SDR for Engineers (PDF) might really be where you want to start reading. Your confusion about what the components of a receiver are is cleared up in the first 15 pages!
Make it through the first four chapters, and you can draw better diagrams than I did here, and that you've got so far. It's literature time that will definitely pay off, because this sounds like you've got a project to approach – and that will require you to make informed decisions :) And: it's actually not all that hard.

¹ and especially not the cases you should be learning first


Changing the cutoff of an IIR lowpass filter with minimal transients?

How about running two filters parallel with each other. Both filters get the same input signal but only one filter will output the filtered signal because there is a crossfader behind them which will pass only one of the filters output.
You can change the parameters of the currently offline filter, and fade the crossfader to it. Parameter changing without any transient.


Here's a paper with a pretty thorough explanation of proper scaling in order to obtain the power spectrum or power spectral density. Essentially N should be the ENBW or the sum of the window squared.


An ideal LC circuit has no resistance, and hence no resisitive loss mechanism. Therefore if it starts to oscillate (due to, for example, a non-zero initial voltage or current stored in the L or C), then it will oscillate indefinetely on its own, yielding a situation as described in your question.
Zero input response of a linear system is the response under no input excitation but due to non-zero initial conditions.
Note the distinction between the terms linear time-invariant (LTI) and the broader term linear alone. The discussion in the Lathi's book is about linear systems in the more general sense.
Note also that a linear (but not necessarily LTI) system can have non-zero initial conditions, and hence zero-input response; but a strictly LTI system cannot have non-zero initial conditions, because one property of LTI systems is the initial rest property. Hence zero-input response of LTI systems is zero by definition.
So by definition a system with zero-input response cannot be LTI, but it can be linear nevertheless...


what is the relationship between a spectrogram and the uncertainty principle heisenberg?

In Quantum Mechanics, state of system specified by a vector (wave function which is a vector in function space), and you could use different basis to represent this vector (imagine one vector in 2 different coordinate system which lead to different components in those coordinate system but they both represent the same vector). For systems composed from moving particles, one basis is the position basis and the second one is momentum basis. The transform between these two basis is the Fourier transform. Considering the wave function give us a probabilistic interpretation, to find the momentum or position of a particle you have to find the average and after that to have a sense about error, you have to find the standard deviation.
The Heisenberg Uncertainty principle, tell us the multiplication of position error and momentum error could not be smaller than some value, for any possible state of system.
Considering this is a property of Fourier transform you could extend that to signals where the signal amplitude and it's spectrum are different representation of same thing in different basis, and say multiplication of effective bandwidth of pulse (as standard deviation of power spectrum around central frequency) and the effective pulse width (as standard deviation of signal's power around it's center in time) could not be smaller than some value for every possible signal.


Low pass filter with different non linear phases

Design your desired magnitude response
Make it minimum phase
Cascade with any type of allpass filter to create your desired phase response


BIBO stable LTI system frequency response for this input signal?

For a sinusoidal input $x(t)=\sin(\omega_0t+\varphi)$, the output of a real-valued LTI system is given by

$$y(t)=|H(j\omega_0)|\sin(\omega_0t+\varphi+\arg\{H(j\omega_0)\})\tag{1}$$

where $H(j\omega)$ is the system's frequency response. The amplitude of the output is scaled by the magnitude of $H(j\omega)$ evaluated at $\omega=\omega_0$, and the phase is shifted by the argument of $H(j\omega_0)$.

If the input is a constant ($x(t)=c$), then in analogy with $(1)$, the output is given by

$$y(t)=|H(j0)|\cdot c\tag{2}$$

Since the system is linear, we can use superposition, and compute the output for each input component and then sum up all the individual components. So the output corresponding to the input $x(t)=c+\sin(\omega_0t+\varphi)$ is given by

$$y(t)=|H(j0)|\cdot c+|H(j\omega_0)|\sin(\omega_0t+\varphi+\arg\{H(j\omega_0)\})\tag{3}$$

So there is no contradiction between a gain of $6$ at frequency $\omega_0$ and a gain of $0$ at frequency $0$ (DC). From the given input and output you can deduce the values of $H(0)$ and $H(\omega_0)$.


Sampling $x(t)=\cos(4\pi t)+\cos(2\pi t)$

How can we know the used sample frequency fs by looking at y(t)?

You can't.

The sampling theorem states that any sampling higher than twice the highest signal frequency allows for perfect reconstruction. In you case any sampling frequency higher than 4 would result in the same $y(t)$.


What is the correct block diagram of CDMA communication system?

Since spreading and symbol modulation are memoryless linear operations, the order makes absolutely no difference. Both methods of writing the same thing down are equally valid, and the one you choose mainly depends on what is convenient for the rest of your text.


Filter long audio data using cascading biquad filter

You are not storing the filter states between buffers.

In this case, each filter is a 2nd order filter, which has 2 state variables (the values stored in the 2 delays implementing each filter).

So for each new frame, you are starting with state 0, introducing a glitch between frames.

You must store the resulting state for each filter, and apply it for each next frame.

Built-in Function: [Y, SF] = filter (B, A, X, SI)

 If the fourth argument SI is provided, it is taken as the initial
 state of the system and the final state is returned as SF.  The
 state vector is a column vector whose length is equal to the
 length of the longest coefficient vector minus one.  If SI is not
 supplied, the initial state vector is set to all zeros.

So for each filter, you must store the corresponding SF, and provide as SI on the next frame.


Extrapolate a 2D array using Fourier Transform

You can use the Fourier transform for a lot of things – it can be part of a particular step in a scaling algorithm. Still, what you're looking for are scaling algorithms.

In this particular case, your original data size is (74 × 128) – what that means in wavelengths or any physical unit doesn't matter to algorithms – and you want to scale it to (512 × 512).

So what you'd need to get 128 -> 512 would be simply inserting 3 zero rows after every image row and then low-pass filtering the result. A 128-FFT, followed by a zero-padding the 512-128 higher bins and then doing a 512-IFFT would do that – but not like you'd want it to, but with circular effects at the border; what you'd be doing mathematically would be a cyclic convolution with a 1/4 band-wide filter.

Same goes for the 74 -> 512 dimension – but here you'd need to interpolate higher first (to 37·512, to be exact, i.e. by a factor of 37·256) just to decimate down back to 512 again. And then you'd still be stuck with the circular convolution.

What you in fact want to do is apply resampling by a rational rate – the core filter for that could be applied by fast convolution, i.e. with the help of overlapping multiple FFTs, but it's not really the point of the algorithm (just like "using a screwdriver" isn't the central part in "building an airplane").

So, look for any of the existing image scaling algorithms.


The spectral efficiency depends on the pulse shape. The basepand BPSK signal can be written as $$s(t)=\sum_k a_k p(t-kT_b),$$ where $a_k$ is equal to either $\sqrt{E_b}$ or $-\sqrt{E_b}$, $E_b$ is the bit energy, $T_b$ is the bit interval (so that the bit rate is $R_b=1/T_b$), and $p(t)$ is a Nyquist pulse. The bandwidth of $s(t)$ is equal to the bandwidth of $p(t)$.

Unfortunately, many textbooks insist on limiting the discussion to square pulses, where $$p(t)=\begin{cases}1,\text{ if $0\leq t <T_b$}\\0,\text{ otherwise.}\end{cases}$$ This means that the spectrum of $s(t)$ is infinite, strictly speaking. This would seem to preclude any discussion of spectral efficiency.

However, when bandwidth is understood as an engineering (not mathematical) concept, we see that we can truncate the spectrum to a desired bandwidth $B$ and, if the pulse $p(t)$ is not distorted too much, communication is still possible.

Different authors make different choices for $B$. As you've seen, some say $B=2R_b$ and others say $B=R_b$. Personally, I use $B=5R_b$ (I'm not as optimistic as most book authors). Fortunately, you can make up your mind: using Matlab (or similar), you can create a close-to-ideal BPSK signal $s(t)$ by sampling it much faster than required, and then filter it using different values of $B$. You can see the resulting distortion and intersymbol interference, and decide what value of $B$ you're comfortable with.

If you have the time and inclination, you can do the experiment with an analog signal generator and analog low-pass filters; it's more work but it feels more "real" than a Matlab simulation.


finding the output of these three systems when connected in series?

You are right and the manual is wrong.

Given $S_1, S_2, S_3$ and the respective input-output signals as below :

$$\\x[n] \rightarrow \boxed{S_1} \rightarrow v[n] \rightarrow \boxed{S_2} \rightarrow w[n] \rightarrow \boxed{S_3}\rightarrow y[n] \\$$

$\\$

$ x[n] \xrightarrow{S1} v[n] = \begin{cases}
x[n/2] & \text{if n is even;}\\
0 & \text{if n is odd.}\\
\end{cases}$

$ v[n] \xrightarrow{S2} w[n] = v[n] + (1/2)v[n-1] + (1/4)v[n-2] $

$ w[n] \xrightarrow{S3} y[n] = w[2n] \\ \\$

You would compute the output $y[n]$ in three steps as: (you have compined step 1 into 2 actually)

$v[n] \xrightarrow{S2} w[n] $

$
\begin{align}
w[n] &= v[n] + (1/2)v[n-1] + (1/4)v[n-2]  \\ \\
&=  \begin{cases} {
x[n/2] + (1/4)x[(n-2)/2] ~~~~~~ \text{if n is even;}\\
(1/2)x[(n-1)/2] ~~~~~~~~~~~~~~~~~~~~~~ \text {if n is odd.} }
\end{cases}
\end{align}
$

$\\$

and then, feeding $w[n]$ to $S_3$, we have:

$w[n] \xrightarrow{S3} y[n] = w[2n] \\$

$ y[n] =  \begin{cases}
x[2n/2] + (1/4)x[(2n-2)/2] & \text{if 2n is even;}\\
(1/2)x[(2n-1)/2] & \text {if 2n is odd.}\\
\end{cases}$

Hence (your) the solution is right:
$$y[n] = x[n] + (1/4)x[n-1]. $$

solution manual is wrong:
$$x[n] \xrightarrow{S1,S2,S3} y[n] = x[n] + (1/2)x[n-1] + (1/4)x[n-2].$$

You can verify this simply by putting some test signals, or even easier by writing a simple MATLAB/Octave/Python script to implement the $S_1,S_2,S_3$ system.


How to calculate the upsampling rate (sample per symbol) for a recorded raw IQ .wav file?

It is a unit conversion. You have the sampling rate given, $F_s$, and it has units of $\frac{\text{samples}}{\text{second}}$. You also have the symbol rate, $F_{sym}$, and it has units of $\frac{\text{symbols}}{\text{second}}$. You are trying to figure out what is the $\frac{\text{samples}}{\text{symbol}}$ so it is this division:
\begin{align}
\frac{\text{samples}}{\text{symbol}} &= \frac{ F_s }{ F_{sym} } \\
&= \frac{ \big(\frac{\text{samples}}{\text{second}} \big)}{ \big(\frac{\text{symbols}}{\text{second}}
\big) } \\
&= \frac{\text{samples}}{\text{symbol}} \checkmark
\end{align}


Does MSAA make sense in audio / 1d data streams?

No. The reason is that MSAA works in the visual domain. Optically, we have hardwired edge detectors. MSAA eliminates unintended edges in generated images.

Acoustically, we do not have such edge detectors, and therefore there's no reason to remove them. Also, loudspeakers already filter out such edges mechanically.


The sinogram does make sense if the object is deformable, although if it was sagging under the force of gravity, it should be registering maximum deviation when viewed at 90 degrees (gantry 0 degees is "North"). This one seems to be registering maximum deviation at approximately 135 degrees.

If the glass tube can bend then it can compress as well. If we forget true orientation for a bit and focus on the top sinogram, then the top part of the tube is compressing while the bottom part is "expanding" to create the bend.

In the bottom image, you seem to be correcting in a "naive" way. That is by shifting columns up or down (vertically) by some amount dictated by the curvature of the walls of the tube.

However, the object is deformable. Which means that the slice that is vertical to the gantry's plane of rotation is NOT vertical to the plane of deformation. In other words, shifting up or down distorts the image because it picks up adjacent "voxels" that do not belong to the slice you are trying to correct for. Another way to look at this is that you are decreasing the resolution at the compressed areas and increasing it at the expanded areas.

This last point is important when applying corrections because it will insert some error that is proportional to the elasticity of the tube / specimen. That is, geometrically, the objects might reconstruct at more accurate locations but their density "measurements" might be slightly "blurred".

If you absolutely have to correct the sinogram then you would be looking at image warping and specifically polynomial (or spline) warping.

The approach is to fit a deformable rectangular fine grid over the deformed object and then REMAP the deformed grid on to a parallel grid using interpolation (to fill in those gaps because of curvature I am hinting above). In other words, you will be scanning the (deformed) sinogram, left-right, top-bottom ALONG the curve of deformation and each one of these points will be REMAPPED on to the straight grid locations. (Something similar to this)

This assumes that the contents of the tube are homogeneous (which is probably not true). If they are not, then their deformation depends on their elasticity, the tube deforms differently than the specimen and while you will be correcting one, you will be distorting the other.

Anyway, what you can try to do is to fit a polynomial on the top boundary (corresponding to angle) and at each point use vertical (to the slope of the polynomial) segments (corresponding to y) to recover the warping polynomials.

For more information regarding the background of polynomial image warping please see this link (Sections 2, 2.6) and this link. For more information on using MATLAB to recover and apply these transforms please see this link and more specifically this link (or possibly this link too).

Finally, an idea would be to set the tube upright and rescan but that orientation has its own set of distortions and compromises.

Hope this helps.


MATLAB: How to copy column vectors from a matrix into a cell array?


The frequency response of a CIC filter is identical to a moving average filter:  The response approaches a Sinc function as the number of samples in the moving average increases. (And likewise for higher order CIC filters, $CIC^n$ will be identical to cascading n moving average filters). You mentioned decimating which will cause the response I mentioned to fold into your final Nyquist band according to the decimation chosen.
More importantly, you mentioned observing the DAC output. The DAC output will also have a Sinc frequency response with the first null at the sampling rate (for a simple non-jnterpolating DAC), so not sure if you are getting confused by seeing that? This is because the DAC is a zero-order hold of your desired digital samples (staircase output prior to any filtering). You can view that process as a convolution of your desired samples with a pulse that has a width equal to your sampling rate $F_s$; thus you will have a multiplication in Frequency of your desired digital spectrum with a Sinc with the first null at $1/F_s$.

Below are some figures I have showing this that may offer more insight. A CIC can be used as a low pass filter, or in a multi-rate system as a decimating filter or interpolating filter. The filter is convenient for multi-rate since it is very simple in structure and will have nulls at all the frequency centers where aliasing (decimation) and images (interpolation) occurs, at the expense of a "passband" droop that can be easily compensated for.  (For details on a simple 3-tap CIC compensator, see how to make CIC compensation filter)

The plot below shows basic CIC (Cascade Integrator-Comb) filters without including the optional decimation or interpolation.  This is handy to see as the impulse response of all the structures shown are IDENTICAL. The CIC is a reformed moving average filter (so you can get the frequency response easily in Matlab using freqz[1 1 1 1] for example for a 4 point moving average such as that shown below, or compute it yourself using $1+z^{-1}+z^{-2}+z^{-3}$, where $z=e^{j\omega}$ with $\omega$ going from 0 to $2\pi$ (the unit circle).  That is the frequency response of a CIC (before considering decimation). As the number of samples increase in the moving average, this response approaches a Sinc function. Anything less and we get a Sinc fucntion with all the frequencies above Nyquist folded in (as one way to view it).  The feedback structure is an accumulator, which is a digital integrator, and the longer delay with a subtraction is a "Comb" filter, since on its own it's frequency response in frequency will cycle up and down mimicking the appearance of a comb.  As shown the two structures can be applied in any order (if not changing rate with decimating or interpolating) and will match the moving average filter.

To see how the CIC filter is identical to a moving average filter, you can compare both with a series of impulses (unit samples) at the input which is very insightful, or directly from the math below.

$$1+z^{-1}+z^{-2}+z^{-3} + ...z^{-(n-1)} = \frac{1}{(1-z^{-1})}(1-z^{-n})$$

For example for our case in the figure above:

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

which is:

$$(1-z^{-1})(1+z^{-1}+z^{-2}+z^{-3}) = (1-z^{-4})$$

And then it is fairly simple to multiply out the polynomial and prove equivalence.

Also observing the poles and zeros of the filter structure is insightful: A moving average filter will have n-1 zeros symmetrically spaced around the unit circle with the omission of a zero at z=1 (DC) (but spaced with that position), where n is the number of taps in the filter, coinciding with the Sinc frequency response we expect. For example, the 4 tap filter shown will have zeros at j, -1 and -j.  The comb filter using a delay of n samples will have n zeros symmetrically spaced around the unit circle including the z=1 position, so for the filter shown with n=4, this would be 0, j, -1, -j. The accumulator has a pole at z=1. So cascading the accumulator (integrator) with the comb results in the same zeros as the moving average filter since the zero at z=1 gets cancelled with the pole!

Finally to show the beauty of the CIC, I will add decimation (interpolation is similar).  Observe the simplification that occurs when we combine a decimator with a comb structure. Since we are throwing away three out of every four samples (in the case of the decimate by 4 shown below), we can reverse the order of comb and decimation as shown and still get the identical result.  (Mathematically this is directly derived using the Noble Identities).

So combining these relationships results in the CIC Decimator. For any decimation structure a anti-alias filter is required prior to decimation (for the same reason that an anti-alias filter is required prior to A/D conversion-- a decimator is simply the all digital version of an A/D converter-- it is a digital-digital converter and will create aliasing just as an A/D converter does, for the same reasons- consider the analog signal as a digital signal with an infinite sampling rate if that makes it any easier to see).

What is so great about the CIC is that it is completely reconfigurable just by changing the decimation rate. The filter response of the anti-alias filter will automatically adapt to be the moving average filter with a length equal to the decimation rate, which conveniently places zeros at the center of all the alias locations!  This is demonstrated in my figure I pasted below. This was from CDMA hence the frequency axis with 630KHz passband and "ChipX" shown but demonstrates well the operation of a decimate by 4: all the images shown are frequencies that will fold into the desired passband if not rejected prior to decimation.

Also here is an interesting history of the CIC filter that Rick Lyons had written: https://www.dsprelated.com/showarticle/160.php


PortAudio: how to generate multiple sounds

The beauty of PCM (pulse-coded modulation, i.e. time samples of audio) is that the mathematical sum of two signals sounds like the acoustical superposition.

So, just add your two sines sample-wise.


Demonstrate that $x(t)$ and $\hat{x}(t)$ are orthogonal

You have to use the following form of Parseval's identity:

$$\int_{-\infty}^{\infty}x(t)y^*(t)dt=\int_{-\infty}^{\infty}X(f)Y^*(f)df\tag{1}$$

where $^*$ denotes complex conjugation, and where $X(f)$ and $Y(f)$ are the Fourier transforms of $x(t)$ and $y(t)$, respectively. With $y(t)=\hat{x}(t)$ and $Y(f)=-j\,\textrm{sgn}(f)X(f)$, we get from $(1)$

$$\int_{-\infty}^{\infty}x(t)\hat{x}^*(t)dt=j\int_{-\infty}^{\infty}|X(f)|^2\textrm{sgn}(f)df\tag{2}$$

For real-valued $x(t)$, $|X(f)|^2$ is an even function, and, consequently, the integrand on the right-hand side of $(2)$ is odd. This implies that the integral on the right-hand side of $(2)$ equals zero if it is interpreted as a Cauchy principal value.

Note that we've used the fact that $x(t)$ is real-valued. For complex-valued $x(t)$, the integrals in $(2)$ are generally not zero.


How to calculate the power of a finite length signal?

If the signal $r(t)$ is of finite length (and if integral of its square during its definition interval is also finite) then its average power given by
$$ \bar{P_x} = \lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2} |r(t)|^2 dt = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{L} |r(t)|^2 dt = \lim_{T \to \infty} \frac{[\text{constant}]}{T} = 0$$ will be zero...

Btw, you may compute an instantaneous power of the signal, simply as the square of it $|r(t)|^2$...


What type of sensors have complex Gaussian noise?

No physical system suffers from, or produces, complex Gaussian noise. Quadrature systems can be modelled as being complex, with the in-phase branch corresponding to a real value, and the quadrature branch corresponding to an imaginary value. If each branch presents Gaussian (thermal) noise, then the noise can be modelled as complex Gaussian, too. However, this is just a modelling technique: the actual physical noise is never complex.

To answer your first question: sensors that operate in quadrature, and are modelled using the complex envelope, can be modelled as having complex noise.


A signal containing no pure sinusoidal components

It contains no spikes in its frequency spectrum. (What is the Fourier transform of a sinusoid?)


How do I simulate a purely real signal being received on multiple antennae in an array?

The input at each array will be complex, because the wireless channel "spreads" the energy of the in-phase and quadrature parts of the signal into each other.

One simple way to think about this is transmitting a pure in-phase carrier: $$s(t) = \cos(2\pi f_c t).$$ The wireless channel will add a random phase, so the received signal is $$r(t) = \cos(2\pi f_c t + \phi)$$ (ignoring changes to the signal amplitude for the moment). Note that, since $\phi$ is random and changes over time, there is no real way for the receiver to determine it.

Now calculate the signal out of the quadrature branch of the receiver: $$ r_Q(t) = \text{LPF} \lbrace r(t) \sin(2\pi f_c t + \theta) \rbrace,$$ where $\theta$ is the phase of the receiver oscillator. If you do the math, you'll see that the energy of $r_Q(t)$ is not zero.


Derive Frequency Representation of Impulse Train Function

@MattL suggested a nice, simple way to see the above result.

But if you want to see the result in the normal analysis equations you mentioned, you can do like below.

Say S(t) is a periodic train of impulses.So S(t) can be written as

$$\ S(t)= \sum_{n=-\infty}^{\infty} \delta(t-nT)$$

Now if you take the fourier series of S(t),you can write S(t) as

$$S(t) =\sum_{n=-\infty}^{\infty} C_ke^{jnw_ot} $$

Where $C_n$ are exponential fourier series coefficients and $w_o$ is the fundamental frequency .

So from exponential fourier series we know that

$$C_n= (1/T)\int_{-T/2}^{T/2} S(t)e^{-jnw_ot} dt$$

Now in the above expression substitute the value of S(t) from the first expression.

So $$C_n = (1/T)\sum_{n=-\infty}^{\infty}\int_{-T/2}^{T/2} \delta(t-nT)e^{-jnw_ot} dt$$

Now, you have to make an observation, if you observe the integral, it's from -T/2 to +T/2. During this integral period, observe that only a single impulse $\delta(t)$ exists.All the other impulse functions in the summation occur after T/2 or before -T/2. So in total the above equation for $C_n$ can be written as

$$C_n = (1/T)\sum_{n=-\infty}^{\infty} \delta(t)e^{-jnw_ot}$$

From sifting property we can write the above as

$$C_n = (1/T)e^{-jw_on(0)} = (1/T)$$

Now put this value of $C_n$ in the first S(t) equation

$$S(t) = (1/T)\sum_{n=-\infty}^{\infty} e^{jnw_ot}$$

Now find the fourier transform of above equation

$$1\Longleftrightarrow 2\pi\delta(w)$$

$$e^{jw_ot}\Longleftrightarrow 2\pi\delta(w-w_o)$$

So the fourier transform is
$$S(jw) = (2\pi/T)\sum_{n=-\infty}^{\infty} \delta(w-nw_o)$$

This should help.


Digital Lowpass Filters, basic questions

If you want as little delay as possible, I recommend notch filters at your resonance frequencies. If your resonance frequencies are fixed, it should be easy. Otherwise, a low-pass filter at about 100 Hz should do the trick.
No your delay depends mostly on the cut-off frequency, the order of the filter and the filter type. An order-2 IIR filter will have a smaller delay than a 30-tap FIR filter . Your sample rate will not affect (or barely) the delay.
Typically, the lower the cut-off frequency, the higher the delay and vice versa. Also, an order-2 low-pass IIR will typically have a smaller delay than an order-4 low-pass IIR for the same cut-off frequency.
If you can live with non-constant group delay, an IIR should do the trick.
Yes it does, but it assumes that you know what your resonance frequency is. You can maybe adapt your notch filter to the resonance frequency in real-time.. it is feasible but can be a lot of work.


Train Cascade Object Detector example in MATLAB

You are not doing anything wrong.  This particular example is there for a quick illustration, and it does not produce a very good detector. It trains a 5-stage detector from a very small training set.  A decent detector may have 20 stages or more, and you would need thousands of positive samples and negative images to train it.

Take a look at this tutorial for more information.


Typically the goal of crest factor minimization is to find the time signal with the small-ish crest factor for a given amplitude spectrum in the frequency domain. If your goal is different, please update your question.

The only "dial" you have is the phase spectrum. Hence the primary output of a CF minimization process is a phase spectrum. With the exception of a few special cases there are no good analytical solutions for this (at least as far as I know), so most methods are iterative or heuristic.

If you choose a random phase, you typically end up with a crest factor of about 4 to 5, which is simply a consequence of the central limit theorem. Your amplitude distribution will be roughly Gaussian, which in theory has infinite crest factor. However, in practice, amplitudes higher than 5 are extremely unlikely (around 1e-6) so they rarely occur or the error resulting by simply clipping them off of negligible

Iterative methods work the following way:

Initialize with a random phase
Go to time domain
Clip time domain signal slightly below the target crest factor
Go back to the frequency domain and repair the amplitude spectrum
Repeat until you reached your target crest factor or you time out

The idea here, is that every iteration step drives the phase spectrum closer to one which represents a lower crest factor. This works quite well with reasonably "broad" spectra and you may be able to get to about 1.2 or so.

Sparsely populated spectra are more tricky since there are less degrees of freedom.

An alternative method is to construct the time domain signal directly by concatenating chunks of sine waves where the dwell time at each frequency is proportional to the target energy at this frequency. Then repair any spectral errors in the frequency domain. This will create more "sweep like" signals and tends to works better for spectrally sparse signals. CFs will come about 1.6 or so. In the special case of a white signal this is equivalent to the "squared phase" approach


No it's not correct.

Derive it like this :

1-) $2N$-point DFT of $x[n]$ is: $X_{2N}[k] = \sum_{n=0}^{2N-1} x[n] e^{-j \frac{2\pi}{2N} k n }$

2-) $Y[k] = X_{2N}[2k+1]$ be the odd-indexed samples of $X_{2N}[k]$

3-) We are looking for $y[n] = \text{N-point IDFT}\{ Y[k]\}$.

4-) Elaborate on step-2 and step-1 to see that $Y[k] = \text{N-point DFT}\{ x[n] e^{-j \frac{\pi}{N}n}\}$

5-) Then from steps 3 and 4 we get :

$$y[n] = \text{N-point IDFT}\{ \text{N-point DFT}\{ x[n] e^{-j \frac{\pi}{N}n}\} \}  $$

$$ y[n] = x[n] e^{-j \frac{\pi}{N}n}  $$


Algorithms for finding fundamental frequency based on ACF result

im assuming your signals are audio and harmonic and you are not doing multiple F0 estimation where you have multiple tones played simultaneously. in that case there are a few easier methods. the normalisation of short time autocorrelation function. r(lag) = (1/N)*sum(x(n)*x(n+lag)
N is length of x(n)
F0 here is obtained from the inverse of the lag that corresponds to the maximum of r(lag) within a predefined range (for example your range could be from 20 to 250, hence at lag 196 you have a maximum, again play around with this range). to avoid detecting an integer multiple of the period, short lags are better than long ones.
here is a C function for this, probably worth playing around with.

void autocorr(
    float *r,       /* (output) autocorrelation vector */
    const float *x, /* (input) data vector */
    int N,          /* length of data vector */
    int order       /* largest lag for calculated
                      autocorrelations */
)
{
    int     lag, n;
    float   sum;

    for (lag = 0; lag <= order; lag++) {
        sum = 0;
        for (n = 0; n < N - lag; n++) {
            sum += x[n] * x[n+lag];
        }
        r[lag] = sum/N;
    }
}

An alternative method for detecting peaks, is center clipping. after your autocorr function you could code something like:

y(n) = a(n) - C,  if a(n) >= C

y(n) = a(n) + C, if a(n) <= -C

y(n) = 0 otherwise.

C is some percentage of the max amplitude of the autocorr vector.

all i have showed you above is the very basic stuff, since F0 analysis can be very complex (probably why people have not answered). you can find more in speech codec books or music transcription books.


Calculating Stability of Phaser Feedback Loop

A feedback loop like this (with enough phase diversity) will be stable if the magnitude of the loop gain is less than unity. In this specific example, it's really simple: make sure your feedback coefficient needs to be less than 1.

So it's $|z|<1$, x and y don't really matter.


Once an audio signal has been digitalized, it has no units and it is treated as a number.

You could use dBFS (decibels relative to full scale) as unit, though. The value of 0 dBFS represents the maximum digital signal level.


Are time series data always contain noise?

Periodically taken ticker data is very likely not a set of samples of a band-limited function, thus can contain aliasing artifacts which can be considered noise.  This is possibly a much larger contribution than quantization and time jitter noise.

Absence of noise requires the existence of well defined signal being measured.  Proving that such a signal really exists is a separate and likely far more difficult problem for some types of time series.


Downsampling - how many extra samples to include

I want to downsample a finite section of an equidistant time-series
  by doing an FFT and dropping the upper part of the spectrum.

Why? This one of the least efficient way to down sampling.  See Why is it a bad idea to filter by zeroing out FFT bins?

In order to down-sample you need to low pass filter first. There are many different ways of building a low pass filter and almost all of them are better than zero-ing FFT bins.

  What heuristics could lead me to choose a reasonably appropriate number of extra samples (probably by specifying an accuracy parameter)?

You do this by first specifying the relevant requirements for your application: signal to noise ratio, phase distortion, pass band ripple, stop band ripple, latency, CPU/memory constraints. Some of those may be very important, some of those may be irrelevant, but it's different for each application.  Then you design a filter that makes sure you meet your requirements.


When did signal processing become a field of its own?

There's a nice historical account given by the IEEE Signal Processing Society, that references 1948 as the year the field began. I'm not sure if that's necessarily accepted as "official" by other sources, but their reasoning is pretty convincing (at least as a ballpark):

  The year 1948 may be regarded as the annus mirabilis in the emergence
  of the discipline of signal processing. For in that year Claude
  Shannon published the epoch-making "A mathematical theory of
  communication"; Bernard Oliver, John Pierce, and Claude Shannon
  published the classic argument for the use of pulse-code modulation;
  modern digital methods of spectrum estimation were introduced;
  error-correcting codes were introduced; audio engineering achieved a
  new prominence; and the IEEE Signal Processing Society, albeit under a
  different name, was established.

That quote is from the beginning of the second chapter, Annus Mirabilis: 1948.


$$X_1(z) = \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}} =\frac{A}{1 - 1.5z^{-1}} + \frac{B}{1 - 0.8z^{-1}} \tag{1}$$

  Mistake is that $A \neq 2$ and $B \neq -0.5$.

Well then, what are they ? There aren't any $A,B$ where the above equation is true. Why? Take the limit $z \rightarrow \infty$. The term $\frac{A}{1 - 1.5z^{-1}} + \frac{B}{1 - 0.8z^{-1}}$ factors to maximally $z^{-1}$ in the numerator whereas $\frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}} $ contains a $z^{-2}$ term!!!!. So a contradiction !! The correct decomposition is as follows:
$$X_1(z) = \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}} =\frac{A}{1 - 1.5z^{-1}} + \frac{B+Cz^{-1}}{1 - 0.8z^{-1}}$$

Taking a common denominator in $\frac{A}{1 - 1.5z^{-1}} + \frac{B}{1 - 0.8z^{-1}}$, equation $(1)$ is now equivalent to
$$2.6z^{-1} - 1.8z^{-2} =A+B + (C -1.5B - 0.8)z^{-1} -1.5Cz^{-1}$$
By identification, we get the following system of equations
$$\begin{bmatrix}
    1.0000 &   1.0000   &      0\\
         0 &  -1.5000  &  1.0000\\
         0  &       0  & -1.5000
\end{bmatrix}
\begin{bmatrix}
A\\B\\C
\end{bmatrix}
=
\begin{bmatrix}
         0\\
    3.4000\\
   -1.8000
\end{bmatrix}
$$
which will give you
\begin{align}
A &= 1.4667\\
B &=  -1.4667 \\
C &= 1.2
\end{align}
Can you continue ?


Help with audio EQ cookbook BPF filters and Q

i believe that what you're seeing is the compression of bandwidth done by the frequency warping inherent with the bilinear transform.  now, if you look at cookbook, the bandwidth $BW$ is expressed in terms of octaves.  so, for your spec:

$$ BW = \log_2(20250) - \log_2(19750) = 0.036069 \text{ octave} $$

try using that number for $BW$ or the corresponding $Q$ as shown in the cookbook:

$$ \frac{1}{Q} = 2 \sinh\left( \frac{\ln(2)}{2} BW \frac{\omega_0}{\sin(\omega_0)} \right) $$

i come up with $Q = 7.63359$.

your formula relating bandwidth in linear frequency to Q is accurate for analog filters but is accurate for digital filters (designed using bilinear transform) only for frequencies much lower than Nyquist.  the reason why is that frequency warping of the BLT moves frequency points downward and even after compensating the center frequency, it compresses the bandwidth in comparison to the analog counterpart.  and 20 kHz is quite close to Nyquist, so we expect this frequency warping and bandwidth compression.

note that for analog, the equation relating $Q$ and $BW$ (in octaves) is slightly different:

$$ \frac{1}{Q} = 2 \sinh\left( \frac{\ln(2)}{2} BW \right) $$

note that the factor $\frac{\omega_0}{\sin(\omega_0)}$ increases the effect of bandwidth, so the digital bandwidth need be smaller than the analog bandwidth especially as $\omega_0 \to \pi = $ Nyquist.


How to compute the energy of two signals?

You simply have to apply the definition of energy. Assuming all signals involved are real, the energy of $g_3(t)$ is given by

\begin{align}
E &= \int_{-\infty}^\infty g_3^2(t) \, dt \\
&= \int_{-\infty}^\infty \big(g_1(t) + g_2(t)\big)^2 \, dt \\
&= \int_{-\infty}^\infty g_1^2(t) \, dt + 2\int_{-\infty}^\infty g_1(t)g_2(t) \, dt + \int_{-\infty}^\infty g_2^2(t) \, dt.
\end{align}

One interesting case is where $g_1(t)$ and $g_2(t)$ are orthogonal, in which case the middle integral is zero. In that case, the energy of $g_3(t)$ is the energy of $g_1(t)$ plus the energy of $g_2(t)$. One simple example is two rectangular pulses that do not overlap.


How do I use the step and ramp functions to create a specific signal?

You have a confilct in your definition for t < 0.

I'm going to assume you mean f(t) = 0 for t < -1, with the rest of the definition unchanged.

The simplest way I can think of doing this is to use u(t) as a tool to switch functions on and off over the desired ranges of time.  For example, if you want a function to be on over a range of time, say t1 to t2, multply the function by u(t-t1)-u(t-t2).  This turns the function on at t1 and off again at t2.

One representation of your function is:

$$f(t) = [u(t+1)-u(t-1)]t + 2u(t-1)$$

For this representation, we've used two functions: t and the constant function 2.
The function t is switch on at t1=-1 and off at t2=1.  The constant function '2' is switched on at t2=1 and never switched off.

If you need to use the ramp function instead of t, your function is:

$$f(t) = -[u(t+1)-u(t)]r(-t) + [u(t)-u(t-1)]r(t) + 2u(t-1)$$

Here we've used three functions: -r(-t), r(t) and 2.


My first and second choices are Python and IDL, not necessarily in any order.

Python is popular with scientists, engineers and artists, well-known, on all platforms, free and open.  You need Numpy and PIL (Python Imaging Library).  Python is more and more commonly used in many places where Matlab used to rule.

IDL from IIT-VIS (I think it's a different name now) is similar to Matlab and the like, but nicer syntax for multidimensional arrays. Very popular with astronomers, geologists, and various others.  It has flaws.  IDL costs $$ but luckily there's a free version named Fawlty Language. It's not open source, but is free (last time I checked.) I don't have the link handy, but google it.   There's also GDL which is a Gnu attempt to create an IDL-compatible interpreter, but last time I tried it, it wasn't so good.

I don't like Octave because it's clumsy for three dimensional and higher arrays.  Images are only 2D, but I often need to process a stack of them, and also may have RGB or other types of channels.


Which pattern (circle pattern or checkerboard pattern) should be used for automotive camera calibration (fisheye/wide webcam)?

The answer here suggests that checkerboard patterns may yield more accurate (subpixel) calibration results and be more robust.

You may have edited your question because the title asks which pattern but the text asks about rows and columns. In either case, you may consider using checkerboard instead of circle pattern maybe?


Is There a Sparse Representation for Noise?

Let's think about it in a different way - Generate Noise from a Dictionary.

Let's create a Dictionary $ A \in \mathbb{R}^{m \times n} $ where each of its rows is normalized (Has Euclidean Norm of $ 1 $) and generated by a Gaussian Random Vector.

Now, let's create $ N $ random vector $ {\left\{ {r}_{i} \right\}}_{i = 1}^{N} $ by:

$$ {r}_{i} = A {g}_{i} $$

Where $ {g}_{i} $ are random vectors $ {g}_{i} \in \mathbb{R}^{n} $ with only $ k \ll n $ non zero elements (Let's say they are pre defined and deterministic otherwise we'll have to deal with multiplication of random variables)..

Now, the output vectors are valid Gaussian Random Variables and given enough of them and some properties of $ A $ you'd be able to build a Dictionary which they are sparse with relation to it.

Yet it means the dictionary is random by itself and next time when you get new set of random data it won't be able to represent it correctly (It doesn't generalize).

While the above is not rigorous (It means to develop some intuition, not analytic analysis) it tells you why when we deal with finite amount of data it can be done but for noise as in general it won't work.


Characterization of transfer functions with no local peaks

As far as I know there is no simple characterization of filters with a monotonic magnitude response. I could only come up with a weak time domain characterization of LTI systems with monotonically decreasing amplitude for positive frequencies (cf. Eq. $(6)$ below).

Let $h(t)$ be the real-valued impulse response of an LTI system, and let $H(\omega)=A(\omega)e^{j\phi(\omega)}$ be its Fourier transform. Furthermore, define $A(\omega)$ and $\phi(\omega)$ such that $A(\omega)\ge 0$ is satisfied. Now we require

$$\frac{dA(\omega)}{d\omega}\le 0,\qquad\omega>0\tag{1}$$

Let $r_h(t)$ be the deterministic auto-correlation function of $h(t)$:

$$r_h(t)=\int_{-\infty}^{\infty}h(\tau)h(\tau+t)d\tau\tag{2}$$

Note that the Fourier transform of $r_h(t)$ is given by $A^2(\omega)$. By a basic property of the Fourier transform we have

$$\begin{align}t\,r_h(t)&=\mathcal{F}^{-1}\left\{j\frac{dA^2(\omega)}{d\omega}\right\}\\&=\frac{j}{2\pi}\int_{-\infty}^{\infty}\frac{dA^2(\omega)}{d\omega}e^{j\omega t}d\omega\\&=-\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{dA^2(\omega)}{d\omega}\sin(\omega t)d\omega\\&=-\frac{1}{\pi}\int_{0}^{\infty}\frac{dA^2(\omega)}{d\omega}\sin(\omega t)d\omega\tag{3}\end{align}$$

Note that since $A(\omega)\ge 0$, condition $(1)$ implies

$$\frac{dA^2(\omega)}{d\omega}\le 0,\qquad\omega>0\tag{4}$$

From $(3)$ and $(4)$ we get the following inequality:

$$|t\, r_h(t)|\le -\frac{1}{\pi}\int_{0}^{\infty}\frac{dA^2(\omega)}{d\omega}d\omega=\frac{A^2(0)-A^2(\infty)}{\pi}\tag{5}$$

which finally results in

$$|r_h(t)|\le \frac{A^2(0)-A^2(\infty)}{\pi |t|},\qquad t\neq 0\tag{6}$$

In sum, condition $(1)$ implies that the deterministic auto-correlation of the corresponding impulse response must satisfy inequality $(6)$.


How to Linearly Combine Two Unbiased Estimators of One Parameter without Knowledge of Their CoVariance?

Basically I'd do something like:

$$ \frac{ {\sigma}_{2}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2}} {r}_{1} + \frac{ {\sigma}_{1}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2}} {r}_{2} $$

This is the optimal weighing given knowledge of the Variance only (Well, linear).
Basically it assumes the cross correlation is 0.

Derivation

There are many ways to derive this.
The easiest would be intuition of knowing that Variance is the Linear Measure of the quality of an estimator. But we're leaving intuition our for this moment.

So easiest formal derivation could be using the ideas @Mark Leeds talked about.

Our model is:

$$ {r}_{3} = {w}_{1} {r}_{1} + {w}_{2} {r}_{2} $$

We want to guarantee the mean of $ {r}_{3} $ is un biased so:

$$ \mathbb{E} \left[ {r}_{3} \right] = \mathbb{E} \left[ {r}_{1} \right] = \mathbb{E} \left[ {r}_{2} \right] = {w}_{1} + \mathbb{E} \left[ {r}_{1} \right] + {w}_{2} + \mathbb{E} \left[ {r}_{2} \right] \Rightarrow {w}_{1} + {w}_{2} = 1 $$

So we'll use $ {w}_{2} = 1 - {w}_{1} $. Now we want to minimize its variance. Since it is not biased we're after minimizing its second moment:

$$\begin{align*}
\mathbb{E} \left[ {r}_{3}^{2} \right] & = \mathbb{E} \left[ { \left( {w}_{1} {r}_{1} + \left( 1 - {w}_{1} \right) {r}_{2} \right) }^{2} \right] \\
& = {w}_{1}^{2} \mathbb{E} \left[ {r}_{1}^{2} \right] + { \left( 1 - {w}_{1} \right) }^{2} \mathbb{E} \left[ {r}_{2}^{2} \right] + 2 {w}_{1} { \left( 1 - {w}_{1} \right) } \mathbb{E} \left[ {r}_{1} {r}_{2} \right] \\
& = {w}_{1}^{2} {\sigma}_{1}^{2} + {\left( 1- {w}_{1} \right)}^{2} {\sigma}_{2}^{2}
\end{align*}$$

We assumed $ \mathbb{E} \left[ {r}_{1} {r}_{2} \right] = 0 $.

Now, this is a quadratic form so to minimize just derive it and find what vanishes it:

$$ {w}_{1} {\sigma}_{1}^{2} - {\sigma}_{2}^{2} + {w}_{1} {\sigma}_{2}^{2} = 0 \Rightarrow {w}_{1} = \frac{ {\sigma}_{2}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2} }, \; {w}_{2} = \frac{ {\sigma}_{1}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2} } $$

Pay attention that the solution doesn't require us to know the exact variances of the estimators being fused. Only their ratio.

So in the case above where $ {\sigma}_{1}^{2} : {\sigma}_{2}^{2} = 0.8 : 0.2 $ the solution is $ {w}_{1} = 0.2 $ and $ {w}_{2} = 0.8 $. This is the best you can do without the knowledge of the co variance.

This fused veritable will have variance $ {\sigma}_{3}^{2} = {0.2}^{2} {\sigma}_{1}^{2} + {0.8}^{2} {\sigma}_{2}^{2} = 0.16 \hat{r} $

Wikipedia Solution

You may, as a reference, look at Wikipedia - Sensor Fusion.
They wrote the same equation in a real clumsy manner:

Yet they are equivalent (Just multiply their result by $ \frac{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2}} $)

Pay attention that the sensor model assumes that 2 sensors are not correlated (Well, it assumes noise is property of the measurement equipment, hence this assumption is made).


“Subtracting” one audio stream from another audio stream

This setup shares some similarities with system identification problems, where $A$ would be the input of the LTI system you want to estimate the transfer function of, $B$ being the output; and the "ambient sound" being the additive noise. The LTI assumption is reasonable provided your converters/amplifiers/transducers are of decent quality.

So the steps would be:

Use a system identification technique to find the FIR filter $\hat{h}$ that minimizes the mean-square error between $\hat{h} \star A$ and $B$. A simple method, which might not be the most suitable here, is to divide the cross-correlation of $A$ and $B$ by the autocorrelation of $A$. Explanation here. The limitation is that it will not work well for long recordings (you might better compute your estimates on shorter segments and average them) - and that music is not the best "probe" signal to send into a system to estimate its response.
You can now use $\hat{h} \star A$ as an estimate of the original signal A as "heard" by the microphone and subtract it from $B$ to retrieve the ambient sound.

I gave a shot at this using a music clip (A), applying a reverb and slight amp model to simulate a speaker in a room, then mixing in a cat audio sample to get (B), then estimating an impulse response from the (A, B) pair, then subtracting the filtered A from B. This shows some results but a better FIR estimation technique might help here! (note that I truncated the estimated IR to its first 5000 samples to speed up computations).

Note that there are algorithms for doing this adaptively (such as LMS). This might be more suitable for your problem if $A$ and $B$ are processed in realtime rather than offline. Such algorithms form the basis of echo cancellation systems used in telecommunications.


OFDM Modulation: what's the output of the FFT?

not missing anything. You're right, if you have an N-point (I)FFT, the input is an N-element vector, and the output is an N-element vector.
The wikipedia diagram just "hides" the fact that these N output elements are "emitted" one after the other from the IFFT block.


Basically I am working on the compressed response(poin target) of an image. Here in image i have a point target(white bright dot) and in zoom your can see its RGB cells. I want to compute 2 Dimensional SNR of point target. Point target is the signal and rest of it is noise.
This processed image has an array of (512x1024)rows*coulums= 524288 pixels
My approach is follow the basic formula SNR= 10log(Psignal/Pnoise)
The point target (white bright point) has a main lobe(compressed response) (bump)
The rest are side lobes.
##Power_signal= Main_loob Energy/pixels covered by point target##Power_noise= Rest side lobes/(Total Pixels - pixels of point target)##
I am working on this approach and not completed yet...
I need guidance on 2D SNR:

Is my approach correct
2)Do we consider back ground as noise or not
3)Any other methods that will make the problem easy to solve


Computing the decimation ratio between two m-sequences

I finally asked this question on mathoverflow.net and got a very helpful answer.

The short answer is to follow essentially the following procedure:

Find any root $y_1 \in F_1$ of $p_2(y)$ in the polynomial ring $F_1[y]$ --- these are the polynomials $p(y)$ which have coefficients that are elements in $F_1$ and are themselves polynomials.
Take the discrete logarithm to determine $j$ such that $x^j = y_1$
Any decimation ratio in the cyclotomic coset of $j$ is a valid decimation ratio.

I've written up a longer description in my blog post which includes the use of a factoring method that partitions $p(y)$ as $p(y) = \gcd(p(y),S(x^ky))\gcd(p(y),S(x^ky)-1)$ where $S(u) = u + u^2 + u^4 + u^8 + \ldots + u^{2^N-1}$, which is mentioned in Lidl + Niederreiter and also in http://blog.fkraiem.org/2013/12/01/polynomial-factorisation-over-finite-fields-part-3-final-splitting-cantor-zassenhaus-in-odd-characteristic/.


At each pixel $p$, this $2\times 2$ matrix summarizes the local behavior  of the image (structure) inside a neighborhood of pixel $p$ (a neighborhood weighted by function $w$.

From this matrix, different quantities can be computed, that provide interesting attributes for image analysis. Structure tensors, or second-order matrices, have been used for image segmentation, especially with geodesic or level-set methods, as they can provide insights to image textures or orientations.


What exactly is Savitzky-Golay differentiation filter?

Your first suggestion is correct, the derivatives of the local polynomials are being sampled. From the horse's mouth (Savitzky & Golay 1964):

  Again, if we restrict ourselves to evaluating the function only at the
  center point of a set of equally spaced observations, then there esist
  sets of convoluting integers for the first derivative as well. (These
  actually evaluate the derivative of the least squares best function.)

About your second question: No, differentiation cannot be moved to the discrete-time domain without affecting the result. The number of taps is equal between any given Savitzky–Golay filter and all Savitzky–Golay derivative filters for the same polynomial degree and number of input samples used to construct the polynomial. This indicates that Savitzky–Golay derivative filters are not derived from the corresponding Savitzky–Golay filter by further discrete-time filtering, because that would increase the length of the impulse response.

As a reality check, it can be observed that the zeros of the frequency responses of the 0th and 1st derivative filters do not match, so one can not be produced by filtering the other:

Figure 1. 2nd-order 9-point Savitzky–Golay (blue) and 1st derivative filter (red) magnitude frequency responses.

Octave source code for the plot:

graphics_toolkit("gnuplot");
pkg load signal;
N=9;
f = sgolay(2, N, 0)((N+1)/2,:);
g = sgolay(2, N, 1)((N+1)/2,:);
[hf, wf] = freqz(f, [1], 4096);
[hg, wg] = freqz(g, [1], 4096);
plot([wf, wg], 20*log(abs([hf, hg]))/log(10));
xlim([0, pi]);
ylim([-100, 0]);
xlabel('\omega');
ylabel('abs(frequency response) (dB)');

References

Abraham. Savitzky and M. J. E. Golay, "Smoothing and Differentiation of Data by Simplified Least Squares Procedures", Analytical Chemistry 1964 36 (8), 1627-1639
DOI: 10.1021/ac60214a047


What is the inverse Z transform of this:

You just need to multiply the terms and sort by powers of $z$. If you do so you should get

$$X(z)=\frac{z^2}{4}-\frac54+z^{-2}\tag{1}$$

So the work is in multiplying out the terms, not in actually figuring out the inverse $\mathcal{Z}$-transform, because from $(1)$ you can immediately write down $x[n]$:

$$x[n]=\frac14\delta[n+2]-\frac54\delta[n]+\delta[n-2]\tag{2}$$


Alternatives to offline Kalman filtering

If you used a commercial GPS receiver to obtain your tracking data, the raw trilateration data from satellites are already processed "using Kalman filtering (more common) or least-squares (less common) estimation algorithms" (see https://insidegnss-com.exactdn.com/wp-content/uploads/2018/01/marapr13-Solutions.pdf). Pay attention also to the counterintuitive statement that the article proves in the section *No Silver Bullet": the differences between least-squares and Kalman are surprisingly minor for the generic application in GPS receivers.
With this in mind, additional post-processing is quite unlikely to improve your data set, but instead your data can deteriorate as a result of this post-processing. It seems you cannot denoise your data any better, but you can analyze your route (whether the highway route passes "canyons", sharp corners or the like) and, knowing the features of your GPS receiver taken from the manufacturer's datasheet, estimate the precision of your tracking data.
The article referenced also gives a partial answer to your highlighted question about an "alternative".


Direct sequence spread spectrum modulation

Multiplying $r_t$ by the pseudo random spreading code does different things to $m_t$ and $i_t$: Because $m_t$ is correlated to the spreading code, it is despread while keeping its power, so the power spectral density will be higher in a small bandwidth (due to remaining modulation). On the other hand, $i_t$ is not correlated and is translated to wideband noise again.

This gives you the opportunity to filter with a bandpass and just use the band, where the signal to noise ratio is above unity.

The noise is not eliminated, but its overall power is reduced, because your bandpass rejects some of the noise (that is still wideband) and passes all of the signal (that is now narrowband).


How to remove salt & pepper noise from speech signal?

If you are forced to use the Discrete Haar Transform, then note that the Discrete Haar Mother Wavelet is nothing but a [1 -1], or in other words, a discrete differentiator, as shown below:

As such, it can/will severely accentuate sudden transitions in your waveform, ie, salt and pepper noise.

Note, that the first dyadic scale of the DHT will be composed of coefficients computed from the signal correlated with the Haar Mother wavelet. At this level, salt&pepper noise will correlate very well, and will give you large wavelet co-efficients in the wavelet space. One can then simply null the large values out, and then perform an inverse Discrete Haar Wavelet Transform, (IDHT).

It is strange that you would be required to use the Discrete Haar Wavelet Transform, (over, say, the Stationary Wavelet Transform), because nulling co-efficients in the DWT space and inverse transforming can lead to aliasing. This effect might be small for the problem you describe however.

Another possibility would be to simply perform the DWT, identify the large wavelet co-efficients within the last half of the DWT transform, (corresponding to high correlations between SP noise and the Haar Mother Wavelet), identify those large co-efficients, and simply null them out in the time domain, without doing an inversion. This will get around the problem of aliasing in the inverse transform.


Digital Anti Aliasing Filter For Waveform Band Limiting

I am generating my waves in a "Raw" mathematical way, meaning that i am creating a ramp for a saw wave and my square wave consists of pure 1's and 0's. The issue is that this technique is inherently not band limited. I would like to implement low pass filter into the algorithm in order to filter of the harmonic being generated that are higher than the nyquist limit.

You can't do this.  You've already aliased the frequencies when you generated the waves in that way.  You can't remove them with a filter.

You're essentially sampling the pure mathematical signal without using an antialiasing filter first.  It's not possible to add the antialiasing filter after sampling.  You will just filter the already-aliased signal.  Bandlimiting has to occur before sampling.  You need to change the way you generate the waveforms if you want them to be bandlimited.  You could generate them using additive synthesis, as I've demonstrated here, or you could generate them the simple way at a much higher sampling rate to minimize the aliasing, then antialias filter, then decimate.  Or you could use BLIT or BLEP, etc.  (I think those are both "reduced-aliasing" not 100% eliminated, but I'm not sure.)


Mathematically speaking, is a "signal" a function or the set of outputs from a function

A signal is a physical quantity (e.g. voltage) carrying information, or a set of values (e.g. samples in discrete case) of the given function for different values of the underlying independent variable.

In the diagram above $\epsilon(t)$ is a mathematical modeling (i.e. function) that characterizes the set of inputs (i.e. scalers) going into the controller.

Strictly speaking, a signal is the "observed" (measurable) physical quantity while a function is a mathematical model (abstract representation) describing the relationship between the involved variables.


Don't understand coefficients in MATLAB "butter" low pass result

You are using wrong syntax. You should call freqz function as:

freqz(b,a)

You do not divide any coefficients. In case of non-recursive filters (FIR) you use only coefficients stored in b vector. When you are using the recursive LCCDE (Linear Constant Coefficients Difference Equation) (mostly IIR filters) then you need also the a vector. I suggest you to set a to 1 if you use non-recursive FIR filters.

You should also watch out for numerical problems with transfer function representation as mentioned in final paragraph of this manual entry and this answer. It's always safer to use z-p-k representation.


Playing sounds with only one frequency at a time

Yes, as long as we are also allowed to vary the amplitude versus time and we're allowed any time period in which to change frequency. Any single valued waveform can also be described in terms of instantaneous frequency versus time, where the instantaneous frequency is the derivative of phase so in this case the rate of change in phase at any given point in time.
If the OP is referring to FM modulation itself (encoding the instantaneous amplitude of a real waveform to the frequency of a constant envelope sinusoid), then yes this is also done where the parameter that can be varied is the frequency deviation (linear mapping of waveform amplitude to frequency) resulting in the modulation index and ultimately the bandwidth required of the modulated signal.  This modulated signal would need to be demodulated with a frequency discriminator to recover the original waveform.


Applying Viterbi algorithm to compensate for ISI in PSK31

(After some deliberation, I have a guess at answering my own question. Someone tell me if I'm wrong.)

The first step is to look at the contribution of each pulse at the decision point. The filter is normalized so a single pulse results in 1. And in this particular case, the tails of adjacent pulses work out to 1/6.

We can then look at any possible combination of three pulses, summing their contributions together. If we say the three pulses are $x_{n-1},\ x_{n},\ x_{n+1}$, then the output $y_n$ at the decision point is:

$$ y_n = 1/6\:x_{n-1} + x_{n} + 1/6\:x_{n+1}$$

To give some examples:

111: 1/6 + 1 + 1/6 = 1.33
010: -1/6 + 1 + -1/6 = 0.67
110: 1/6 + 1 + -1/6 = 1

Say state 01 means the previous bit was a 0, and the one before a 1. We can then ask, "if we add another bit, what does the output become?" For example, if the next bit is a zero that yields 001, with an output of -1 according to the formulation above.

We can then build a state table:

╔═══════╦═══════╦════════╗
║ input ║ state ║ output ║
╠═══════╬═══════╬════════╣
║     0 ║    00 ║ -1.33  ║
║     1 ║    00 ║ -1     ║
║     0 ║    01 ║ -1     ║
║     1 ║    01 ║ -0.67  ║
║     0 ║    10 ║ 0.67   ║
║     1 ║    10 ║ 1      ║
║     0 ║    11 ║ 1      ║
║     1 ║    11 ║ 1.33   ║
╚═══════╩═══════╩════════╝

And the associated trellis diagram:

The notation "0/-1.33" means "to send a 0 from this state, the output is -1.33".

Either the magnitude or the square of the difference between the observed values and the expected values can be used as the metric. The Viterbi algorithm will find the path through the trellis with the lowest metric regardless of what metric is used, as long as it's something that monotonically increases with the distance between observation and expected.

The Viterbi algorithm can then be followed in the usual way to determine the most likely sequence that explains the received signal, knowing how adjacent pulses will influence the received values.

I've found these videos especially helpful in understanding the evaluation of the Viterbi algorithm:

Viterbi Algorithm, by AppBooke Xperience
Viterbi Algorithm, by Keith Chugg


The difference is the expected input. The cosd function expects the input to be expressed in degrees, and the cos function expects the input to be in radians. So you have cosd(theta)==cos(deg2rad(theta)). This is just from the MATLAB documentation page.

Are you using these correctly? For example, you do cosd(2*pi*F*t). So we have $\frac{2\pi \text{ radians}}{\text{cycle}}.\frac{F \text{ cycles}}{\text{second}}.[\text{seconds}]=2\pi F \text{ radians}$, as the input. You should be using cos not cosd because your input is expressed in radians.

BTW, I do not check the correctness of the actual code. That is on you.


A better explaination for 'there's no flipping of the impulse response'

Flipping the impulse response is really just a matter of perspective.  The LTI system doesn't care about perspective.  In any case, here is a graphic showing a system that takes an input of color weighted impulses.


How to get main frequency from digital signal

Generally, yes, you should perform FFT and take the highest peak from the first half of the transformed signal. However, if the frequency your are looking for is constant (or close to constant), i.e. doesn't change (much) throughout the observation time, I would recommend using Welch's method instead of a simple FFT.
If the accuracy is not good enough you can perform interpolation around the peak, or on the entire signal using zero-padding.


Lowpass Butterworth filter design (own poles' calculations)

Ok, I've solved this problem myself. I've used equation from this article mentioned in first post. I've just implement it in matlab not enough precise.

clc, clear all, close all;

%ANALOG BUTTERWORTH FILTER DESING
%https://engineering.purdue.edu/~ee538/AnalogFilterDesigns.pdf
%Butterworth LP with -3dB at fc and N order and G gain

wc = 1000;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain

%analog filter design
%poles and zeros

%poles searching
for k = 0:N-1
   sk(k+1) = wc*exp(j*(pi/2))*exp(j*(2*k+1)*(pi/(2*N)));
end

[NUM,DEN] = zp2tf([],sk,G);
b0 = DEN(end)*G      %gain normalization

Hs = tf(b0,DEN)

figure();
freqs(b0,DEN)
figure();
plot(sk./10^2, 'x')
xlim([-10,10]);
title('Poles placement in s-plane')

Many thanks for your help,
Maks


time downsampling vs. frequency downsampling

Your second equation suggests that $y[n]$ is obtained by performing an $N$-point inverse DFT of an $M$-point DFT $X[k]$ of an $M$-point sequence $x[n]$.

The result can be interpreted as shifted copies of $x[n]$ are superimposed, causing possible aliasing in time, which depends on whether $ N \geq M$ or $N < M$.

If $N \geq M$, there is no aliasing and $y[n]$ will be an $N$ point sequence whose first $M$ samples will be identical to $x[n]$ for $n=0,1,...,M-1$ and remaining $N-M$ samples will be zero-padded.

$$
y[n] = \begin{cases} { x[n] ~~~,~~~ n = 0,1,...,M-1 \\ ~~~~~0 ~~~,~~~ n= M,M+1,...N-1} \end{cases}
$$

Otherwise if $N < M$, then $y[n]$ will be a time-aliased version of $x[n]$, as shifted copies of $x[n]$ would be overlapping each other. If you want any non-aliased samples in $y[n]$, then $N > M/2$ in general and in such a case $y[n]$ could be defined as:

$$
y[n] = \begin{cases} { x[n]+x[n+N] ~~~, n = 0,1,...,M-N \\ x[n] ~~~~~~~~~~~~~~~~~~~~~,~~~ n= M-N+1,...,N-1 } \end{cases}
$$


A question about the Commutativity Property in discrete signals block diagrams

If you apply the multiplication with $b_0$ at the output, then that multiplier is after the feedback loop. Consequently, if the output (after the multiplier) is $y[n]$ then the input to the delay element is $y[n]/b_0$. The corresponding difference equation is

$$y[n]=b_0\left(x[n]-\frac{a_1}{b_0}y[n-1]\right)=b_0x[n]-a_1y[n-1]\tag{1}$$

which is identical to the one with the multiplier at the input, as expected.


Some thoughts on it at least.

First, if you can shown linearity and time-invariance and you know that the filter has the correct impulse response you are home. Given that the filter is stable. So in this case there is no need to run different other input signals and the frequency response is given based on the impulse response.

Checking impulse response is obviously quite simple. Linearity and time-invariance I guess is actually really complicated (any specific input signal you test may just work out of coincidence, but others may not). However, as the filter code in most cases shouldn't be that complicated, it should be possible to eye-ball the code and see that it does what is expected (no non-linear operators, no if-cases based on signal values etc).

Remains stability. For IIR filters this is a really complicated issue and I will not go into that discussion. For FIR filters, they will always be stable (unless they are implemented using some recursive algorithm). However, you may run into numerical issues relating to overflow (most likely not underflow) and round-off noise. You can find formal proofs that an FIR filter will not overflow (maximum output is the sum of the absolute values of the impulse response coefficients times the maximum input magnitude). For round-off errors you can find statistical expressions, but they are not of much use in formally proving anything (but rather as an argument that the expected output with "infinite precision" and the actual output still conforms).

Also not that all impulse responses etc should be evaluated using the actual (possibly rounded) coefficients not the (possibly higher precision) ones obtained from the filter design tool. This is because the filter will always realize the transfer function given by the coefficients that it use, no matter how they were derived. This may still not be the impulse response you get at the output though, as there may be round-off issues. An easy way to see that is to use a scaled impulse (i.e. not one) and note that typically not only the magnitude of the frequency-response will change, but also the shape. The filter will still implement the same transfer function, but the round-off noise added to the two different input signals are different.

Hope this gives you a bit more insight.


Methods to quantify randomness (or complexity) in a signal

There are two different concepts:

If you think as your signal as a single random variable $X$ that is emitting values, then what you want is to calculate the Entropy of the random variable

http://en.wikipedia.org/wiki/Entropy_estimation

If you are considering the entire random signal or stochastic process, then you have to estimate the autocorrelation function. The most random signal possible is White noise, in which every sample of the signal is not correlated with any other except itself.

http://en.wikipedia.org/wiki/Autocorrelation


HINT:

Set $f(t)=x(5t)$ and use your rule number 2: $$f(t)\star \delta(t-t_0)=f(t-t_0)=\ldots$$


Do energy of signal and power spectrum (in db/hz) indicate towards same signal property?

On a straight forward approach, energy of a signal given by $E = \int \limits_{-\infty}^{+\infty}|x^2(t)|dt$ is a constant number whereas PSD can be obtained from FFT which is a function of frequency. PSD is the square of the absolute value of FFT.

In either case the value of energy or PSD gives an insight of the amplitude/strength of a signal.


FFT of a signal with 0 added between samples

You can use expansion analysis to deduce what happens when a sequence is zero staffed in between its samples.

Let $x[n]$ be your original length $N$ sequence. Then its expansion by $2$ yields the new sequence $y[n]$ of length $2N$ :

$$ y[n] = \begin{cases} { x[n/2] ~~~,~~~n = 2m \\ ~~~~~ 0 ~~~,~~~ \text{otherwise} } \end{cases} $$

Then the DTFT relation between $Y(e^{j\omega})$ and $X(e^{j \omega})$ will be : $$Y(e^{j\omega}) = X(e^{j 2 \omega})$$

Then the DFT relation between N-point $X[k]$ and 2N-point $Y[k]$ will be

$$Y[k] = Y(e^{j \frac{2\pi}{2N} k}) = X(e^{ j 2 \frac{2\pi}{2N} k}) =  X[k] ~~~,~~~ k = 0,1,...,2N-1$$

In other words, $Y[k]$ will be two copies of $X[k]$


Initial conditions for systems described in state space - LTI or not?

If you look at Chapter 5 of:

https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-011-introduction-to-communication-control-and-signal-processing-spring-2010/readings/MIT6_011S10_notes.pdf

which is entitled "Properties of LTI State Space Models" , equation 5.33 doesn't seem to have a problem with initial conditions, or any other book (I stand corrected, there is one book) that I'm aware of.  Unless Oppenheim was touched with insanity, I'm inclined to accept his characterization that initial conditions don't disqualify a LTI system as "not linear" by his use of the term "zero input linear".

At the beginning of the notes, (and in Oppenheim and Shaefer 3rd edition) a LTI system is given as:

$$
y[n]=\sum_{k=-\infty}^{\infty} x[k] h[ n -k]
$$
 which doesn't require $h[n]$ to be causal or stable.  $x[n]$ doesn't have to satisfy $x[n]=0 \;\text{for} \; n <0$ .

There is emphasis  in the text that one needs to consider the entire history of $x[n]$, not just for $n \ge 0$.

let $$x[n] = \hat{x}[n]+\tilde{x}[n] $$
  where
$$ \hat{x}[n] = \begin{cases} x[n]\; \text{for}\; n<0 \; \text{and}\\
0 \; \text{for}\; n\ge 0 \end{cases}$$

and
$$ \tilde{x}[n] = \begin{cases}0\; \text{for} \; n<0 \; \text{and}\; \\ x[n] \; \text{for}\; n\ge 0 \end{cases}$$
 by linearity.

$$
y[n]=\sum_{k=-\infty}^{-1} \hat{x}[k] h[ n -k]  + \sum_{k=0}^{\infty} \tilde{x}[k] h[ n -k]
$$
if $y[n]$ is causal,
$$
y[n]=\underbrace{\sum_{k=-\infty}^{-1} \hat{x}[k] h[ n -k]}_{\text{zero input linear}}  + \underbrace{\sum_{k=0}^{n} \tilde{x}[k] h[ n -k]}_{\text{zero state linear}}, \quad n \ge 0
$$

The essential point is that initial conditions account for prior input.  Where $n=0$ is referenced for $x[n]$ is arbitrary, which is another manifestation of time invariance.  Initial conditions are not some arbitrary value vexing the system.  if $x[n] = 0 $ for $n <0 $ the initial conditions are zero.

Let's try something else. Let $z[n]=\tilde{x}[n+1]$ (advance by one sample) and with $\tilde{x}[n]$, the system was LTI without controversy.
But now,
$$
y[n]=\underbrace{ z[-1] h[n]}_{\text{zero input linear}}  + \underbrace{\sum_{k=0}^{n} z[k] h[ n -k]}_{\text{zero state linear}}, \quad n \ge 0
$$
and now we have an initial condition. A forward shift of 1 sample would make an LTI system nonlinear?

The logical fallacy at the root of question is to use the definition of zero state linearity and applying it to the zero input case.


You just need to binarize the image then find the hull. In Mathematica it looks like this:

<< ComputationalGeometry`
img = EdgeDetect@Rasterize["managers. In\n more lean\n result of\nss\nkforce and\nt \
points and"]
pts = Position[ImageData@img, 1]
pgp = PlanarGraphPlot[pts, ConvexHull@pts][[1, 2]]
Graphics[{Black, Point@pts, pgp}]

If this is inadequate I would find the convex hull of the text two lines at a time, then fuse the results.


Impulse Response with non-delayed feedback

Take the z-transform of your equation to get $y(z) = x(z) h(z) - y(z) h(z)$ which you can simplify to $$ y(z) \left( 1 + h(z) \right) = x(z) h(z) $$ and therefore get the transfer function
$$ H(z) := \frac{y(z)}{x(z)} = \frac{h(z)}{1+h(z)} $$

Now plug in the z-transform $h(z)$ of your response $h(n)$ and perform an inverse z-transform to get the total impulse response $H(n)$.

An approach that is more basic and does not require z-transforms works for your particular example by solving the recursion equation explicitly. The inner system with the response $h(n) = u(n) \left( \frac{1}{2}\right)^n$ can be written with the recursive difference equation $y'(n) = x'(n) + \frac{1}{2} y'(n-1)$ where $x'$ and $y'$ are the output and input of the subsystem. Then we can write the full system in terms of the following equations:

$$ y(n) = y'(n)$$
$$ x'(n) = x(n) - y(n)$$
$$ y'(n) = x'(n) + \frac{1}{2} y'(n-1) $$

Eliminating $x'$ and $y'$ results in

$$ y(n) = x(n) - y(n) + \frac{1}{2} y(n-1) $$

which simplifies to

$$y(n)=\frac{1}{2}x(n)+\frac{1}{4}y(n-1)$$

This system has the impulse response $$H(n) = \frac{1}{2} u(n) \left(\frac{1}{4}\right)^n$$


Unscented Kalman Filter - Multiple Consecutive Measurement Updates

I will use Wikipedia notations - Kalman Filter.
In most models the state transition model matrix $ F $ depends on the interval parameter $ T $. The same goes for the Process Noise Covarinace Matrix $ Q $. For instance, for the constant velocity model:
$$ F = \begin{bmatrix}
1 & T \\
0 & 1
\end{bmatrix}, \; Q = \begin{bmatrix}
\frac{ {T}^4 }{4} & \frac{ {T}^3 }{2} \\
\frac{ {T}^3 }{2} & {T}^{2}
\end{bmatrix} $$
Now, one simple policy you act by, whihc works perfectly in real life and is highly flexible, is as following:

Always remember the time stamp of the last Kalman step, be it prediction or update.
Given new input which is not a measurement do as following:

Take the input time stamp.
Define $ T = {t}_{current} - {t}_{last} $.
Build $ F $, $ Q $ and all other model matrices based on the given $ T $.
Apply prediction step.
Update the time stamp of the last Kalman step to $ {t}_{current} $.

Given new input which is a measurement do as following:

Take the input time stamp.
Define $ T = {t}_{current} - {t}_{last} $.
Build $ F $, $ Q $ and all other model matrices based on the given $ T $.
Apply prediction step.
Apply update step.
Update the time stamp of the last Kalman step to $ {t}_{current} $.

This works and it is very robust. Yet it has the overhead that you need to generate the correct matrices each step.
Also, since the UKF (Unscented Kalman filter) is different only in its Update step, it holds for it as well (For EKF as well).


The appeal of this image is obviously in the numerous lines, which test the aliasing properties of resizing, denoising, and super-resolution algorithms.

It seems Allen Gersho is the source, according to the Acknowledgement section of Embedded image coding using zerotrees of wavelet coefficients

If you want more specific information you could ask him or one of his students from the early 90s, when this image first seems to have appeared.


Difference between sampling rate and baud rate and their relation to channel capacity

First let's estabish the context.
Let's say you want to transmit a sequence of numbers $\lbrace a_k \rbrace_{k=1}^{N} = a_1, a_2, \ldots$ Using PAM, this sequence is transmitted with the signal $$s(t) = \sum_{k=1}^N a_k p(t-kT_p),$$ where the pulses $p(t-kT_p)$ form an orthonormal set, that is, $$\int_{-\infty}^\infty p(t) p(t) \text{d}t = 1$$ (which implies that the pulse $p(t)$ has unit energy) and $$\int_{-\infty}^\infty p(t-jT_p) p(t-kT_p) \text{d}t = 0$$ whenever $k \neq j$, which allows the receiver to recover the transmitted sequence using a correlator or a matched filter.
Second, let's see what Nyquist and Shannon say.
First of all, the pulses in $s(t)$ are transmitted at a rate $R_p = 1/T_p$. This is called the "pulse rate" or the "symbol rate". You want $T_p$ to be short, to transmit information faster. However, this increases the bandwidth $B$ of $s(t)$. The relationship is $$B = \alpha \frac{R_p}{2} = \alpha \frac{1}{2T_p}.$$
The value of $\alpha$ is the "spectral efficiency" of the pulse $p(t)$. The minimum possible value is $\alpha=1$, when $p(t)$ is a sinc pulse. In this case, $B=R_p/2$ or $R_p=2B$. This is what it means to say that "$2B$ pieces of information can be transmitted per Hz of bandwidth available."
If this was all there was to it, then an infinite amount of information could be transmitted in any amount of time. All one has to do is let the $a_k$ be arbitrary real numbers, with infinite precision. Even if one is constrained to choose from a constellation, such as $M$-QAM, one could choose $M$ to be as large as needed.
In practice, however, noise, distortion, and other imperfections prevent transmission of $a_k$ numbers with infinite precision. The simplest model of a communications system assumes that the received signal $r(t)$ is corrupted with white Gaussian noise with power spectral density $N_0/2$: $$r(t) = s(t) + n(t).$$
Transmission under the effects of noise restrict the choice of symbols $a_k$: they must belong to a finite set whose elements cannot be easily mistaken for one another. For example, if we constrain the symbols to take values either $-0.5$ or $0.5$, and the noise magnitude is very unlikely to be larger than $0.5$, then transmission can be done very reliably; however, we pay a price in two different ways:

Power: since the symbol amplitudes cannot be made too small, a minimum amount of power is required;
Rate: since we're limited in the choice of symbol amplitudes, we can't transmit as fast as we wish.

In terms of a $M$-QAM constellation, the consequence is that we can't increase $M$ arbitrarily; when $M$ is made too large, the system cannot transmit information reliably any more.
This is where Shannon comes in. He was interested in the relationship between the entropy of the information source that produces the transmitted sequences $\lbrace a_k \rbrace$, the bandwidth, and the signal to noise ratio in the received signal. He definitely did not say that the maximum bit rate is twice the bandwidth; rather, the maximum bit rate is given by the bandwidth, the SNR, and some specific properties of the sequence $\lbrace a_k \rbrace$. The study of these properties is known as "coding theory".
Finally: is it a coincidence that a signal with bandwidth $B$ must be sampled at rate $2B$, and that a signal with bandwidth $B$ can transmit symbols at rate $2B$? It is indeed not a coincidence. You want to transmit symbols at rate $R_p$; in other words, you want a signal that has specific amplitudes at times $kT_p$:

Think of these amplitudes as signal samples:

The blue signal is the transmitted signal $s(t)$, and since its samples have rate $R_p$, then it has bandwidth $B=R_p/2$.


The explanation in Wikipedia for Top-hat transform - Properties seems pretty simple, for the bottom-hat, another name for the black top-hat.

The black top-hat returns an image, containing the "objects" or
"elements" that:

are "smaller" than the structuring element,
are darker than their surroundings.


Can someone help explain what this FFT workflow is doing to my signal, and why it works?

The answer from Bird has it right.

The described setup is an implementation of a lock-in amplifier.

Lock in amplifiers act as very narrow band filters for detecting the presence of a known signal that is buried in noise.

If you multiply two signals (both simple sine waves) you get a new signal that contains a frequency that is the sum of the original two signals and a frequency that is the difference of the original two signals.  (1000Hz * 200Hz gives an output of 800Hz mixed with 1200Hz.)  This process is used in radio receivers to bring the received radio signal down to a frequency where it is more convenient to do further processing.

Your LabView program makes use of this effect to detect the amplitude and phase of the signal in comparison to a reference signal.

The difference of the input signal and the reference signal is a signal of 0Hz - that's DC.  You get it out of bin 0 in your FFT results.

If the frequency of the input and the reference are not the same, then the difference lands in another bin and is ignored by your detector.

That also makes sense of removing the DC offset from the sample sets - any real DC in the input would cover generated DC from the heterodyning step.

As to why your device is doing that, it is probably so that it can recover a transmitted signal to tell how much of it is coming back to the receiver, and do it pretty much no matter how weak the received signal is or how much noise is present.

Since you say it also takes the phase of the output, it may be using it as part of measuring a distance.

You don't say what the result is used for, or in what conditions you use it so there's not really much more I can say about the "why."


So from the other post your error is:

$$
ERR = \sum \left | G_1(n) - G_2(n) \right|^2
$$

To calculate this should be quite simple. Do an element wise subtraction to get the difference between the two frequency responses at each point. Then square each element (magnitude squared if you have any imaginary components for some reason). Then sum over all the elements to get the total difference/error.

As mentioned in the other post you probably want to normalise this by the reference signal to make detecting abnormality easier.

Similarly you could just take the absolute difference rather than squaring. Squaring is more conventional though as it emphasises big differences more and is related to the variance.


How do I convert a timeseries to a different frequency band?

Mathematically, shifting the frequency of a signal is pretty easy:

Following @OlliNiemitalo's answer, the 0.003 frequency shift can either be done in

time domain, or
frequency domain.

I recommend doing it in time domain, by multiplying the signal with a complex sinusoid,

$$e^{j 2 \pi f_\text{shift} n},\, n \in\{1,2,\ldots\}$$

that way you can get arbitrary frequency shifts.

Still, you're not right about needing to add up the amplitudes -- what you're doing is probably coherent sampling (i.e., the sampling clock doesn't stop after each 2048 sample chunk and starts with a random delay for the next, but continously runs at a stable frequency), so adding up multiple FFTs without reducing them to their amplitude works -- that's the beauty of the DFT being a linear operation, and the bases of the base transform it is being orthogonal.

Think of it like this:

Assume you add up many ($k$) consecutive FFTs. Because the DFT is linear, it'd be the same as adding up many consecutive 2048-vectors of time samples. What happens then is that things that repeat every 2048 samples get multiplied by $k$, whereas things that don't get attenuated or cancel themselves out. That sounds like a bad thing™, however, you have to realize one thing about the frequencies that the DFT can represent in each bin:

Every bin in the 2048-FFT represents one frequency that is an integer multiple of the fundamental frequency (ie. $\frac{f_\text{sample}}{N_\text{FFT}}$). So things that perfectly fall into each bin already inherently repeat every $N_\text{FFT} = 2048$. Hence, doing that will usually (assuming sufficiently stable sampling clock!) simply improve your SNR by $k$.

Concluding, just add up 2048-vectors of your time signal, and transform the result. You will get a 2048-DFT with an increased SNR, which definitely is meaningful.


Any idea what's the OCR tool behind `myfonts.com/WhatTheFont` font recognizer?

I believe this question is off-topic, but this page suggests that the algorithm is described in:


Sending a real number over a noisy channel with limited energy budget

The terms of the kind $N(0, 1)*\alpha/E$ are simply representing noise to the system and the way one defines the multiplicative factor $\alpha/E$ simply changes the variance of the gaussian $N(0,1)$ without impacting the mean.

Note: the mean of $R(E)$ is $X$

So the distribution of $R(E) = X + N(0,1)\alpha/E$ is $N(X,\frac{\alpha^2}{E^2})$ and for the received $R(E) = X + N(0,1)\alpha ln(E)$ is $N(X,\alpha^2ln(E)^2)$.

In the first case as the energy of the signal is increasing the effect of noise is decreasing (it's variance). So it makes more sense compared to second one where the variance on noise is increasing with increasing signal energy. However it does not reflect the true SNR

Edit: The correct way with the update in your question is $\sqrt(E)$ term, this is because this truly represnts the actual SNR, and this is what should be used to decide the minimum error detection rule


How to calculate SNR of dead-band quantizer?

Time align the quantized signal with your input signal and subtract to create an error signal (match the amplitudes such as to minimize the error signal). The SNR can then be determined as the ratio of the standard deviation of the input signal (after amplitude matching) to the standard deviation of the error signal (20Log10 of this ratio to be in dB).

The top diagram that I have shows this for a complex signal example, but it is also equally applicable to real signals.


How to Find the Frequency Response of a Communication Channel from Input and Output Symbols in MATLAB?

I'm not sure about the symbols.
Yet if you have a function in time as the input of the channel and another function in time which is the output of the channel and assuming the channel doesn't have zeros in the spectrum and the input function as well, the channel Transfer Function is given by:

$$ \frac{{X}_{Out}(f)}{{X}_{In}(f)} $$

In MATLAB, just apply FFT on each, and do element by element division.
This will give you the DFT of the transfer function.
Apply IFFT to see the time response function.


Is it possible to decompose an image into basis images using MATLAB?

It doesn't make any difference what image you choose. The basis images are the same. When decomposing a given image into basis images we find the coefficients only and define the given image as a weighted sum of all those basis images.

In the following, you can see the basis images based on 2D DFT of size $8$ ($64$ basis images in total)

Note the different 2-D frequencies. Any given $8\times8$ image is a weighted sum of these 64 images.


I think it means the "apparent" frequency of oscillation. It's poorly worded.

What's happening is that the frequency is increased to the Nyquist rate in the first half and then above it, causing aliasing, in the second.

What looks like the frequency increases in the first half and decreases in the second because of aliasing.


In general case, no, you cannot approximate color of an object from digital image with good accuracy.

Human vision system

What you see is light reflected from the surface of the object. First the object needs to be illuminated with light. This light has some sort of spectrum (probability density function of photon energy). When the light hits the object, the object absorbs some of the light. How much each wavelength of light is absorbed depends on material of the object. Light which is not absorbed, is reflected. Due to the absorption of certain wavelengths, this reflected light has different spectrum than the one which illuminated the object.

You sense the reflected light with four types of light sensitive cells. Each type of cell has its own responsiveness (note the image) for different wavelengths of light. Finally brain combines the data from the light sensitive cells and produces information which you observe as color.

What is color?

So what is the natural value of color? It depends on the reflection and absorption properties of the object, as well as the illumination of the scene, but you probably could describe it as the distribution of observed wavelengths of the light reflected from the object in certain illumination. So if we could reconstruct the wavelength distribution of the observed light, we could describe the natural value of the light color.

How digital cameras capture color?

In digital cameras light sensitive cells are commonly replaced with filter patterns which usually try to capture red, green and blue wavelengths of light. Each filter+light sensitive sensor has certain responsiveness for different wavelengths of light. The responsiveness can vary greatly between different sensor, but they usually try to mimic responsiveness of human vision system (for example see MT9P031 quantum efficiency curves on page 39 of this document).

The RGB image processed by the computer has usually three integer values, red, green and blue. These are weighted integrals of the light spectrum captured by the filtered sensors of the camera.

How to approximate the light spectrum from digital image?

You cannot reconstruct the spectrum from three integrals as there is infinite number of possibilities of light spectrum which would have produced the integrals when weighted with the quantum efficiency of the camera systems.


$x(t)$ is a model for an analog signal. After sampling, the real $t$ variable is turned into an integer variable $k$, corresponding to a regular splitting of the real time interval, starting from some origin time $t_0$, evenly spaced by a sampling period $T_s$ being the inverse of the sampling frequency: $T_s = \tfrac{1}{F_s}$.

A discrete signal is thus:

$$ x[k] = x(t_0+k.T_s) = 2\sin(100\pi (t_0+k.T_s)) = 2\sin(100\pi t_0+k.100\pi/50)) $$

Thus, since $$2\sin(100\pi t_0+k.100\pi/50)) = 2\sin(100\pi t_0+2k\pi))= 2\sin(100\pi t_0))$$

the discretized signal is determined by a constant value, depending only on the initial sampling time $t_0$. You can plot if with many plotting software (even Excel).

The total number of samples is ambiguous. Since $k$ take any integer value, the number of discrete samples is infinite. Since it is constant, $x[n]$ takes only one value.

Explanation: you start from a continuous sine, and only get a constant signal. This is an illustration of the Nyquist-Shannon sampling theorem (or aliasing). If the sampling frequency is not consistent with signal maximal frequency, the discrete signal might be a poor version of the analog signal.


Why am I seeing a negative frequency on FFT graph in Gnu Radio Companion?

The stock FFT (in GNU radio?) is a complex-to-complex transform.  Thus any positive frequency peak you see represents a complex signal (phasor) that can include both real and imaginary components.

Since your cosine (or sine) waveform is strictly real (has zero or no imaginary component), the complex FFT result also includes a negative complex conjugated mirror image.  Since complex conjugation inverts the imaginary component, the two imaginary components (of the positive and negative frequency peaks) cancel out, and thus represent your original strictly real signal (with zero imaginary content).  That is why you see a negative frequency peak.

An asymmetric result can occur if you feed the FFT with a complex IQ signal that has non-zero imaginary components (for instance, baseband SSB).  Thus, if you want to see only a positive frequency peak, then you need to feed your FFT with a positive rotating IQ signal or phasor (cos(t) + i sin(t), IIRC).