Sampling with an alternating phase offset Suppose, I sample a signal $z=\mathrm{e}^{\mathrm{i}\omega t}$ alternatingly in $I=\mathrm{Re}(z(t))$ and $Q=\mathrm{Re}(\mathrm{e}^{\mathrm{i}\alpha}z(t))$. Thus my data array contains
$$\texttt{data =} [I(t_0),Q(t_1),I(t_2),Q(t_3)] \quad \mathrm{with} \quad t_n = n \Delta t = n/f_s$$
For $\alpha=\pi/2$, the Q component is kind of a quadrature component. If I do a Fourier transform of that array, I get 4 peaks in the double sided FFT, one peak at $\omega$ and one at $\frac{f_s}{2}-\omega$ (and in addition the signals at negative frequencies).
If $\alpha$ is not $\pi/2$, I get peaks of different strength.

For $\alpha=0$, I get only two peaks at $\pm \omega$, for $\alpha=\pi$, I get only the peaks at $\pm(\frac{f_s}{2}-\omega)$.
Unfortunately, $\alpha$ is unknown in my case and can take any value. However, it is constant. Is there any good operation that reduces the Fourier transform output to that output that I would get if I would sample with $\alpha=0$ for all time points?
EDIT: After being pointed out in the comments below, I have to admit that the term "quadrature components" might be misleading in this case. It kind of resembles properties of these definitions, although it is not a "true" quadrature component splitting.

Bottom Line
Given a sampled tone at normalized radian frequency $\omega_1$ and real amplitude $A$:
$$x[n] = A e^{j\omega_1n} \label{1} \tag{1}$$
If every other sample is shifted by a constant phase $\phi$, the result in the Fourier Transform will be two tones $x_1$ and $x_2$ (precisely so as the only non-zero values when the sampling rate is an integer multiple of $\omega_1$, otherwise $x_1$ and $x_2$ refer to the maximum bin in vicinity of the two peaks):
$$x_1 = |x_1|e^{j\phi_1}$$
$$x_2 = |x_2|e^{j\phi_2}$$
And $\phi$ is determined using the magnitude and phase of each of the two bins as follows:
$$\phi = 2\sin(\phi_1-\phi_2) \text{atan2}(|x_2|, |x_1|)$$

Further Details
The first tone will be at the bin associated with the same normalized radian frequency $\omega_1$ as the original but shifted in phase by $\phi/2$ and an additional tone will be given due to the phase modulation at the bin associated with the normalized radian frequency $\pi + \omega_1$ and initial phase such that the new tone will always be in quadrature phase with $x[n]$ as $-j (-1)^n x[n]$. The amplitude of the original tone will be reduced according to $A\cos(\phi/2)$ while the amplitude of the additional tone will be $A\sin(\phi/2)$. Given the modified original tone as $x_1[n]$, and the new tone as $x_2[n]$ this relationship to $x[n]$ ahd $\phi$ is:
$$x_1[n] = x[n]\cos(\phi/2)e^{j\phi/2}= A\cos(\phi/2) e^{j(\omega_1n + \phi/2)}\label{2} \tag{2}$$
$$x_2[n] = -j (-1)^nx[n]\sin(\phi/2)e^{j\phi/2} = jA\sin(\phi/2)e^{j((\pi+\omega_1)n+ \phi/2)}\label{3} \tag{3}$$
So given this and not knowing $\phi$, we can determine $\phi$ from the DFT results as follows: Given $x_1$ and $x_2$ as the complex values for two bins where $x_1$ is the frequency bin for the expected tone and $x_2$ is the frequency bin for the interference tone due to the introduced distortion of an unknown fixed phase shift every other sample, with:
$$x_1 = |x_1|e^{j\phi_1}$$
$$x_2 = |x_2|e^{j\phi_2}$$
$x_1$ and $x_2$ will always be in quadrature, and the phase between the $x_1$ and $x_1+x_2$ will always be $\phi/2$. Further, when $x[n]$ has no phase offset as we had introduced in \ref{1}, then the sum  of $x_1+x_2$ will also always be real and the phase can be determined from twice the phase angle of $x_1$ alone. However, what is developed here is the more general case applicable to any starting phase offset for $x[n]$. This general case can be solved for using the following:
$$\phi = 2 \times \text{atan2}(|x_2|\sin(\phi_1-\phi_2), |x_1|)$$
Where $\text{atan2}$ is the 2-element arctangent. The term $\sin(\phi_1-\phi_2)$ will always be either +1 or -1 depending if $\phi$ is positive or negative, since the relationship of the two tones will always be in positive or negative quadrature, so the above can also be written as:
$$\phi = 2\sin(\phi_1-\phi_2) \text{atan2}(|x_2|, |x_1|)$$
To extend this to a real waveform such as in the case of the OP, the original tone from the positive frequency axis would be used as $x_1$, and it's negative image would be used as $x_2$ as depicted in the graphic below:

Using a real waveform that is exactly at normalized frequency $\omega_1 = \pi/2$ would not be resolvable since the image from the added phase distortion and the primary negative frequency  component would both be on the same frequency location.

Derivation:
Relating \ref{1},  \ref{2} and \ref{3}:
$$x[n] = x_1[n] + x_2[n]$$
$$ = x[n]\cos(\phi/2)e^{j\phi/2} -j (-1)^nx[n]\sin(\phi/2)e^{j\phi/2}$$
$$ = x[n]e^{j\phi/2}(\cos(\phi/2) - j(-1)^n\sin(\phi/2))$$
$$ = x[n]e^{j\phi/2}e^{-j\phi/2(-1)^n}$$
$$ = x[n] \text{ for n even}$$
$$ = x[n]e^{j\phi} \text{ for n odd}$$
Intuition:
The above relationships are more intuitive to understand knowing that the Fourier Transform decomposes any single valued time domain analytic function into distinct phasors (as the basis functions) that each must rotate with time with a constant magnitude and at a constant rate. We can see the solution for how this would occur for the OP's problem by doing the decomposition geometrically.
First consider an arbitrary $x[n]$ as a rotating phasor on the complex plane in the time domain as shown on the left, where when sampled the normalized radian frequency corresponds to the phase shift from sample to sample. The Fourier Transform of this on the right is a single impulse in frequency, a positive frequency here given the counter-clockwise rotation of the phasor (positive increasing angle). The single impulse corresponds directly to the single spinning phasor.
x[n] with no distortion:

Now consider if every other sample was shifted forward in phase by a constant amount $\phi$:
x[n] with phase distortion:

To decompose this into phasors of constant magnitude and rate of rotation, the first phasor MUST be the first phasor rotated by $\phi/2$ since the error represents a modulation going back and forth by $\phi/2$ for every sample. And for this same reason the second phasor must be in alternating quadrature relationship such that it is modulating the phase back and forth from sample to sample (resulting in a normalize angular frequency of $\omega_1 - \pi$ when we consider it rotating opposite the original tone, or similarly within the range of $0$ to $2\pi$ at the normalized angular frequency of $\pi+\omega_1$.)
Fourier Decomposition of distorted x[n] into two tones:

And if we look at the relationship between any two samples at a given time instance in the decomposed waveform of two phasors, we see through the geometric relationship how the two must be always in quadrature with the cosine and sine relationship of magnitudes such that the normalized magnitude of the sum of the two phasors will be normalized (to one here) as in the original waveform. Consider the first two samples of the original waveform with the second sample offset by an additional $\phi$ as shown by the red phasors below, together with the first two samples of the decomposed tones showing how the quadrature relationship and relative magnitudes must be maintained:
First Sample of distorted x[n] and the Decomposed Tones:

Second Sample of distorted x[n] and the Decomposed Tones:

In summary, we shift the original tone by $\phi/2$ so that each sample it is midway between the sum of the two decomposed phasors. Then in order for the phase to jump back and forth by exactly $\phi/2$ from sample to sample, the second tone must be in alternate quadrature relationship from sample to sample (so will have rotated $-\pi + \omega_1$ since the original phasor modified with a phase shift of $\phi/2$ still rotates at the same rate as $+\omega_1$). Further for the sum of the two decomposed phasors to be on the unit circle (or whatever magnitude the original tone was), the magnitudes of each are related by the $\cos(\phi/2)$ and $\sin(\phi/2)$ since $\cos(\phi/2) + j\sin(\phi/2) = 1$.
This shows the relationship between the Fourier Transform components in the general case with a single complex tone, and taking the real of the above would result in two symmetric tones with the same relationships as is the case for the OP.  Further, to be precise, the descriptions above are specific to the Discrete Time Fourier Transform (DTFT) in which the time domain waveform is extended for all time and would be a continuous function in frequency with the exact solution shown as 2 weighted impulses in frequency corresponding to the two phasors in the time domain and zero everywhere else. To extend this derivation approach to the DFT, x[n] would be truncated in time (windowed with a rectangular window of duration $T$) which would result in convolving the impulses in frequency with a Dirichlet Kernel (aliased Sinc function) and then sampling the resulting frequency response at a spacing given by $1/T$ where $T$ is the window duration. When the original tone is on a bin center of the DFT (multiple of $1/T$), then the DFT would also result in just two tones for this case. For all other cases, multiple smaller tones would also appear (spectral leakage) but the same relationship would apply between the strongest tone closest to $\omega_1$ and the second tone at $\pi+\omega_1$, where the $N$ frequency bins given as integers $k$ from $0$ to $N-1$ are related to the normalized radian frequency $\omega$ as $\omega = 2 \pi k/ N$ .


Classifying Discrete time filter So let's say i have this frequency response of a digital filter:

My question is how can i classify this type of filter ( low pass, high pass,...)?

It's a bandpass filter whose frequency response suppresses the low and high frequencies (with its zeros at $\omega = 0$ and $\omega = \pi$) and amplifies the band in between.

The bandwidth of the filter can be defined by reference to $-3_{dB}$ points, which would give you approximately $W_{-3dB} \approx 0.2\pi$, with a center frequency of $\omega_c = 0.5\pi$.

I would however not consider this bandwidth to be a descriptive one.


Why is this projection wrong? I am trying to figure out why a 3D cube doesn't project in the right place in the image.

I start with the camera at origin and the cube centered at (0,0,1) world coordinates:

Then I create a 3D rotation R of 45 degrees around the Y axis:

(0.707106781, 0, 0.707106781,
0, 1, 0,
-0.707106781, 0, 0.707106781)

And I set the camera to use rotation R^-1, with a translation of R * (0,0,1) (these are camera to world transformations). The camera should have rotated 45 degrees CCW around the cube (I can confirm I have 0.707106781 translation on x and z).

At this point, I would expect to see the cube with the same scale, still centered, and with one of the edges right in the middle of the image. But this is what I get:

The procedure I use to project a 3d point X is (pseudo code):

bring into camera coordinates: Xc = R^-1 * X - R^1 * translation
transform into image coordinates: x = Xc * focal length / Xc.z + principal point

I have read the theory from various sources and it seems I am following the rules.

This is just artificial data, so not sure if focal length and principal point can have any effect on it? What am I doing wrong? Why is the cube not centered?

EDIT:

I construct the cube like this (duplicate vertices are there because I need them to draw lines):

  cv::Matx31d center(0,0,1);
  double delta = 0.2;

  std::vector<cv::Matx31d> volumeCorners = {
    // back
    (center + cv::Matx31d(-delta, delta, delta)),
    (center + cv::Matx31d(delta, delta, delta)),

    (center + cv::Matx31d(delta, delta, delta)),
    (center + cv::Matx31d(delta, -delta, delta)),

    (center + cv::Matx31d(delta, -delta, delta)),
    (center + cv::Matx31d(-delta, -delta, delta)),

    (center + cv::Matx31d(-delta, -delta, delta)),
    (center + cv::Matx31d(-delta, delta, delta)),

    // right
    (center + cv::Matx31d(delta, delta, delta)),
    (center + cv::Matx31d(delta, delta, -delta)),

    (center + cv::Matx31d(delta, delta, -delta)),
    (center + cv::Matx31d(delta, -delta, -delta)),

    (center + cv::Matx31d(delta, -delta, -delta)),
    (center + cv::Matx31d(delta, -delta, delta)),

    // front
    (center + cv::Matx31d(delta, delta, -delta)),
    (center + cv::Matx31d(-delta, delta, -delta)),

    (center + cv::Matx31d(-delta, -delta, -delta)),
    (center + cv::Matx31d(delta, -delta, -delta)),

    (center + cv::Matx31d(-delta, -delta, -delta)),
    (center + cv::Matx31d(-delta, delta, -delta)),

    // left
    (center + cv::Matx31d(-delta, delta, -delta)),
    (center + cv::Matx31d(-delta, delta, delta)),

    (center + cv::Matx31d(-delta, -delta, delta)),
    (center + cv::Matx31d(-delta, -delta, -delta)),
  };

I use a right-handed coordinate system, with the y axis pointing down, for world and camera coordinates, while image has origin at top left.

My translation should have described a rotation around the center of the cube, but it was wrong. I changed it to:

rotate45 * (translation - center) + center

and I get the cube in the middle with the right angle.


Transformed Direct Form II filter: Z transform: reduce multiplications In the diagram  below there are three multiplication operations required to
compute each filter output sample. I'd like to draw a block diagram of an equivalent filter but requiring fewer than three multiplies per filter output sample.

I need to transform the filter structure to Transformed Direct Form II, that it can be seen that number of multiplication can be reduced to 2. But, how to do that

UPDATE:

As per the filter, are we talking about something like the picture below?
If that's the case, what's up with that unmarked node? What is he, addition? just redirection? I've never seen that before.

I'm trying to extract some generalizable rules about how to construct those neat equations you've made about the output and A & B but I'm not really able to put it into words, because I have only to foggiest notion of what you've done there, would you help me to do that?

It's instructive to look at the actual difference equation implemented by the system. If you define a sequence $w[n]$ as the output of the first (left-hand side) adder you get

$$w[n] = x[n] + Aw[n-1]\\
y[n] = Bw[n] + Bw[n-1]$$

Applying the $\mathcal{Z}$-transform gives

$$W(z) = X(z) + A \ W(z)z^{-1}\\
Y(z) = B \ W(z)\left(1 + z^{-1}\right)$$

which results in

$$W(z) = \frac{X(z)}{1-Az^{-1}}\\
Y(z) = \frac{B(1+z^{-1})}{1-Az^{-1}}X(z)$$

In the time-domain this last equation is equivalent to

$$y[n] = Ay[n-1]+B \left( x[n] + x[n-1] \right) \tag{1}$$

Equation (1) shows that this system can obviously be implemented with only two multiplications. If you want to use the Transposed Direct-Form II structure then just use the structure shown here (bottom figure), set $a_1=-A$, remove the branches with multipliers $a_2$ and $b_2$ and the lower delay element (because you only have a first order filter), remove the multipliers $b_0$ and $b_1$, and multiply the input signal directly with $B$.

EDIT: (in reaction to your updated question)

There are a few mistakes in your figure. Note that in the figure I linked to, the input is on the right-hand side. This might have caused some confusion. Anyway, in your figure where you multiply by $B$, just multiply by $A$ instead. Then remove the multiplication by $-A$ and add a multiplication by $B$ to the input signal $x[n]$, before the path splits. The node you refer to is simply an addition; they obviously forgot the 'plus' sign.

As for the equations, what I wanted to do is to get a difference equation for the output signal $y[n]$ in terms of $x[n]$ and of delayed versions of $x[n]$ and $y[n]$. This shows what the system actually does and in this way the description of the system becomes independent of the chosen filter structure. In order to do that, I used the $\mathcal{Z}$-transform which transforms difference equations into algebraic equations. So in the $\mathcal{Z}$-domain it is straightforward to eliminate the signal $w[n]$ and write $Y(z)$ directly in terms of $X(z)$. Transforming back to the time-domain yields the desired difference equation, which can then be implemented by any structure you like.


Why Is the Canny Edge Detection Used Instead of Sobel / Prewitt Edge Detection Before Hough Transformation? I know that for the Hough Transform to work on an image, it needs to be a binary image. To convert from a grayscale image, an edge detection algorithm should be employed. I notice that people always use Canny edge detection instead of others (Sobel, etc). Why is that?

Canny Edge Detection is considered to be a better (In False Alarm sense) edge detection than those you mentioned.
This is, mainly, due to 2 steps:

Non Maximum Suppression - Edges candidates which are not dominant in their neighborhood aren't considered to be edges.
Hysteresis Process - While moving along the candidates, given a candidate which is in the neighborhood of an edge the threshold is lower.

Those 2 steps reduce the number of "False" edges and hence create a better starting point for farther process like Hough Transformation.


Equalizing frequency response in software I'm playing audio through a speaker which I know to have a non-flat frequency response, with particularly poor response in the low ($<1\textrm{ kHz}$) ranges. My plan to get around this problem is to write a script which processes each audio file before passing it on to the speaker such that the energy is boosted in the ranges where the speaker's response is poor.

I know that it is possible to create a filter which does this (e.g. using the firwin2 function in Scipy), but I can't help but feel that this approach is overly naive, and that I will wind up accidentally introducing unwanted artifacts into the signal when I actually play it through the speaker.

So, is this approach doomed to failure?
If not, then why isn't it something that all speaker manufacturers do already?

Whilst I am not well versed in Scipy, or how the firwin2 function works that you have quoted, this is a very complicated problem that you propose, and there is a lot more to think about than you may realise.

There is one very important component of how sound is radiated by a loudspeaker that is quite often overlooked (especially in marketing specifications for loudspeakers).  Directivity.  You state that your loudspeaker has a particularly poor response below 1 kHz.  I assume that means a “non-linear” response on-axis.  What about off-axis?  This off-axis energy is what interacts with your room, i.e. reflections, and the frequency response in these off-axis directions most likely have a response different to your perfect on-axis response. If you correct just the on-axis loudspeaker response what has that inverse filtering done to the other positions off-axis?  Will you always listen, with one ear, at the intersection of your loudspeaker on-axis directions?

Some loudspeaker manufacturers do actually use technology that attempts to apply “correction” (Genelec, Bose etc…) to overcome room anomalies and positional errors in loudspeaker placement. Whilst the methods employed are very smart, we have to remember that these are just best approximations.  It is also not uncommon to use some form of linear phase FIR filters to correct for minimum-phase anomalies in loudspeaker response for mid-high frequencies, and alternative approaches can be applied for lower frequencies.  Some of these techniques can be better understood by reading text such as: Immersive Audio Signal Processing

Other things to consider: Just how much time data do you window in creating your corrective filter?  And where should you make your measurement? Will you always listen in that position? What I’m emphasising here is that temporal and spatial variation of the soundfield need to be considered also.

So what you are attempting to do isn’t exactly “doomed to failure”, it’s just very complicated, and at the end of the day, the desired result is going to be that which suits your subjective tastes.  I hope that some of the questions I have raised have been helpful, I have raised them simply to remind you that it’s not as simple as inverting a single omnidirectional measurement of a loudspeaker at some point in a room…  Otherwise, yes, every loudspeaker manufacturer would have already done this. Good luck on your endeavour!

I’ll leave a random selection of some references to technical literature related to research in this topic space that you may find helpful:

Combined Quasi-Anechoic and In-Room Equalization of Loudspeaker Responses

Invertibility of a room impulse response

Equalization of loudspeaker response using balanced model truncation

Equalization in an acoustic reverberant environment: robustness results


FFT on a Raspberry Pi I'm trying to determine the processing time required to perform a FFT of length 500,000 on a Raspberry Pi in C/C++. I have had difficulty finding any information that would be useful for determining this. How would I go about doing this?

Since you care about the time it takes, you'll want to use an optimized FFT implementation. That would be FFTw or FFTS, realistically. Historically, FFTw is way dominant in software (I mean, even Matlab uses that), but FFTS works really well on ARM, so might be the better choice for your Pi.
Then you'd just write a minimal C++ program that

creates an fftw_plan  / ffts_plan_t
generates random data (std::uniform_real_distribution is your friend in generating real and imaginary parts of numbers)
uses using sc = std::chrono; auto start = sc::steady_clock::now(); as a "start" timer,
fftw_execute/ffts_execute the plan on your random data
auto finished = sc::steady_clock::now();, and auto ms = sc::duration_cast<sc::duration<double>>(finished-start).count()*1000; to get the execution time in milliseconds.

You might want to run the actual plan multiple times (maybe even modifying the data in between) and divide the total time by the number of executions, to reduce the measurement error.

A note on the raspberry Pi: If you're using a modern RPi, don't use 32 bit operating systems, but go for an Aarch64 OS. In my experience, not having the SIMD instructions available in 64 bit mode really really hurts numeric throughput performance. (For an FFT, the better loop performance of aarch64 compilers might also help significantly; haven't benchmarked that.)


Find reappearing pattern in the signal I am working with a signal that describes energy consumption in some house as a function of time and I am trying to find time intervals that correspond to some particular type of dish washer.

What do we know about this dish washer:

It has 4 heating pulses.
Each pulse has amplitude around 2000 and lasts for at least 120 seconds.
Heating pulses are not periodic.
If there is such a dish washer it was used at least twice in a given period of time.

Question that I am trying to answer:

Was such a dish washer used in this house during given time interval?
If yes, I would like to identify all time intervals when this happened.

Example:

csv with whole signal

Let's say I identified one occurrence of this dish washer, like on this figure and trying to find if there are other similar intervals. At the moment I am doing convolution of the whole signal with this dish washer signal and pick time intervals that give high convolution score. Problem is that signal is additive and if person is using other appliances at the same time I may get high score when there is no dish washer there and low score when it is.

Vladimir: You're trying to detect a certain sample pattern (a certain sequence of sample values) within a longer discrete sequence of input samples. This sounds liked a "matched filter" problem to me. If the pattern of samples you're trying to detect is the samples in the red shaded region of your diagram then performing convolution seems like a smart thing to try first before you explore more exotic signal detection algorithms. By convolution I mean an implementation where you pass your long input sequence of values through a tapped-delay line FIR filter. But be careful! You must make sure that the coefficients of the tapped-delay line filter are a "flipped in time" (reversed in time order) version of the sequence you're trying to detect.


Programming a Vocoder I'm want to program my own vocoder synthesizer like "Songify" but I can't find a simple of vocoder algorithm. Maybe you can explain or tell where to look for information on how vocoder works.

I fear that all answers here are irrelevant to the question.
What is called a vocoder in the music production world has little to do with the phase vocoder used in signal processing. To make matters worse the Songify app referenced by the original post is not an example of vocoder. Let us sort this out!

1. Phase vocoder

The phase vocoder referenced by the other answers is a signal processing technique that can be used to perform time/pitch modification of signals (time-stretching, pitch-shifting), by computing a time-frequency representation of the signal (short term Fourier Transform, or STFT) and then inserting/removing signal frames and then maintaining the coherence of the phase information. Its relationship to voice is only historical, and it is currently used for pitch-shifting and time-stretching in lower-end audio hardware/software. RubberBand is an example of open-source C++ time/pitch alteration library based on a phase-vocoder.

2. Vocoder

When people in the field of music production refer to a Vocoder, they refer to a device extracting the spectral envelope of a signal (usually voice, called the modulator), and filtering an other signal (usually a rich synth texture, called the carrier) with a filter whose response is the extracted spectral envelope. For an example of resulting sound, listen to 0:23 in Kraftwerk Trans Europe Express, or Alan Parsons' Project The Raven from the few first seconds. The resulting effect is a vocal-like timbre applied to the melody or chords played by the carrier signal, giving the feeling that a voice is spoken through a synthesizer.

The vocoder being originally an analog device, it was implemented with two banks of a dozen or more bandpass filters with high Q. The modulator signal is sent through the first filter bank, and the amplitude of all the sub-band signals is tracked with an array of envelope followers. In parallel, the carrier signal is sent through another filter bank. Each sub-band is amplified (with a VCA) with the gains given by the envelope followers. If you read analog, you can have a look at the schematics of a vocoder channel here, from Jurgen Haible's living Vocoder - on top the modulator signal filter, at the bottom the carrier filter and the VCA. Software implementations of vocoders stay close to this, simply because music producers expect vocoders to sound like the classic analog devices! But if you don't want fidelity to the "vintage" devices, and want something cheaper than 40 biquads, another way to achieve the same result is to estimate an all-pole filter (of order 8 to 20 depending on how close you want to get to the original voice) from the modulator signal (AR-modeling) ; and then apply this filter to the carrier. The typical problem here is that you'll need to update your filter coefficients every 20ms frame or so ; so you need a representation of the all-pole filter which handles well abrupt coefficient updates.

3. Auto-tune and pitch-remapping

What Songify does is the following: extract the prosody (pitch contour) of the recorded voice, and alter it so that the resulting pitch contour matches a target melody. This is a bit similar to what auto-tune does, with the difference that auto-tune "rounds" the pitch toward the nearest musically accurate semitone, while Songify just pushes it a target value.

The algorithms at work here are very different from traditional pitch-shifting time-stretching, because the voice signal is monophonic and fits well the source-filter model. Approaches like Time-Domain Pitch-Synchronous-Overlap-Add (TD-PSOLA) are much more efficient, both computationally and in terms of quality, for transparently altering the pitch of voice than generic time-stretch algorithms (usually done with phase-vocoders). These are used for example in speech synthesis to alter the prosody of a synthesized sentence - not unlike Songify indeed! Auto-tune is also based on such time-domain methods (detecting full cycles of the input waveform and resampling them).


Why is modern control theory in the time domain instead of the S-domain? I am not sure if I am missing something, but from what I read I understand that modern control systems work in the time domain rather than the frequency frequency or S-domain.
Wikipedia article snapshot below:

Why is the time domain preferred in modern systems?
I thought the entire point of going to S domain is because it simplifies the equations and associated maths with the Laplace or Fourier transforms.
Is state space representation an alternative to the Laplace transform and is it better/faster than analysis in the S-domain?

I would say that the difference is not about time domain or frequency domain, is more about using transfer function or state space equations.
The state equation for a system with a single input and a single output is
$$\sum_{n=0}^{N} a_n \left(\frac{d}{d t}\right)^n y(t) = \sum_{n=0}^{N} b_n \left(\frac{d}{d t}\right)^n u(t)$$
The state space representation for a system
$$\dot{x}(t) = A \, x(t) + B \, u(t), \\ y(t) = C \, x(t) + D \, u(t)$$
where $y(t)$, $x(t)$ and $u(t)$ all can be vectors, extending this easily to multiple inputs and multiple outputs. For a high order system the number of parameters in state space representation is much larger, and there are multiple possible representations, however there are some canonical representations that are unique to a given system, and may have as few parameters as the transfer function.
Linking the two representations
This can be represented in the Laplace domain as
$$s \cdot X(s) = A \, X(s) + B \, U(s), \\ Y(s) = C \, X(t) + D \, U(s)$$
If you want to understand the dynamic of the system you can manipulate the state update equation as follows
$$s X(s) - A\, X(s) = B \, U(s) \\
 (s \mathbb{I} - A) X(s) = B \, U(s)$$
From the above equation you start to notice that when $X(s)$ is an eigenvector of You have $(s \mathbb{I} - A) v = s v - \lambda v$, and this will define one mode of the system $x(t) = x(0) e^{\lambda t}$, this is why you can study the dynamic of the system in terms of the eigenvalues of $A$.
Proceeding to the transfer function we can solve the above equation for $X(s)$
$X(s) = (s \mathbb{I} - A)^{-1} B \, U(s)$
and replace in the outptut equation
$$Y(s) = C \, X(s) + D \, U(s) \\
C \, (s \mathbb{I} - A)^{-1} B \, U(s) + D \, U(s) \\
(C \, (s \mathbb{I} - A)^{-1} B + D) U(s)
$$
And finally the transfer function is
$$H(s) = \frac{Y(s)}{U(s)} = (C \, (s \mathbb{I} - A)^{-1} B + D) U(s)$$
The result is so difficult to write that I will skip (and that is what is called classical control).
This equation will have rational functions (ratio of polynomials) in $s$, If a system has multiple inputs and multiple outputs it will be a matrix of scalar transfer functions mapping the effect of each input to each output.
Apart from notation compactness, linear algebra methods scale well to higher dimensions. One example of other advantage that I can cite is related to my master's degree. You can study system stability using Lyapunov stability condition, and this leads to linear matrix inequalities that can be applied to uncertain systems (with uncertainty in the model matrices).


Adaptive filter for noise cancellation when measuring some input The context is as follows:  I want to measure (as in, digitize) some input signal $x$ (I have no detailed knowledge about $x$ or its statistical description).

A noise signal $x_N$ contaminates $x$ before I digitize it.  That is, the signal I read is: $x' = x + H_1(x_N)$, where $H_1$ is some unknown filter (let's assume linear, but possibly time-varying).

I have access to $x_N$ (for example, $x_N$ could be the "mains" voltage, or the ripple voltage in my power supply).  I want to cancel the noise component in $x'$ to obtain the best estimate of $x$. This estimate $\hat{x}$ is obtained by subtracting a filtered version of the noise:
$$\hat{x} ~=~ x' - \;H(x_N)$$

where $H$ is my adaptive filter.  The best estimate is that with minimum power (as a function of the filter coefficients), since the noise and the true signal are uncorrelated.

My algorithm aims to minimize $\left(\hat{x}\right)^2$.  At time step $n$, I have:
$$\left(\hat{x}(n)\right)^2 = {\left( x' - \sum_{k=0}^L h_k \, x_N(n-k) \right)}^2$$
To obtain the gradient, I take:
$$\frac{\partial \left(\hat{x}(n)^2\right)}{\partial h_k} = -2 \, \hat{x}(n)\; x_N(n-k)$$
And presumably voilà: each $h_k$ at step $n$ is updated by adding $\mu$ times the above expression to the estimate of $h_k$ at step $n-1$.

Question 1: Does the above make sense?

Question 2: Rather, a concern when I look at the math above --- intuitively, I expect the gradient's magnitude to approach zero as the system gets close to the solution (the optimal filter), but that doesn't seem the case from the equation above: the gradient is a copy of the last $L$ samples of $x_N$ (which in principle has non-zero mean), scaled by the present value of the output (which also, in principle has non-zero mean).

Am I doing something wrong?

I'm a little skeptical of your derivation and notation as it abstains from telling whether a Mean-Square or a Least-Squares metric of error is minimized, but rather just works on an instantaneous error at sample $n$. But I believe you eventually would devise an LMS or RLS filter, hence either of the two would happen.

So, based on your notation the adaptive system error is $\hat{x}(n)$, which is an estimate of the true signal $x(n)$,  and you have two correlated noises $x_N$ and $H(x_N)$ which are uncorrelated with the true signal. The standard noise cancelling setup.

Had you used the steepest descend based quadratic MSE surface minimization procedure, your instantaneous gradient on the error surface would be something like:

$$ \nabla J(w) = - E\{ e(n) \bar{x}[n] \} $$

where the expectation is due to stochastic MSE framework. In your notation this would become:

$$ \nabla J(w) = - E\{ \hat{x}(n) \bar{x}_N[n] \} $$

Now, if the optimum condition is reached, then the error is minimized, and we can assume that $\hat{x}(n)$ converged to the true signal $x(n)$ and you may replace it. Remember this is an approximate conditions actually...

$$ \nabla J(w) = - E\{ x(n) \bar{x}_N[n] \} $$

and since your noise source is uncrorrelated with the true signal, this expectation is zero, yielding the zero gradient at the minimum point.

In your case, you don't specifically use expectations, so might be harder to show why the gradient would go to zero. Furthermore, as always, be aware that gradient would only theoretically go to zero. Any practical implementation will deviate from this depending on how unrealistic your assumptions were.


3-point DFT with two multiplications I'm trying to understand DFT in the context of an integer ring (i.e. the number theoretic transform) and I was wondering how to reduce the number of multiplications by the twiddle factors.
Consider 3-point DFT:
X_0 = x_0 + x_1       + x_2
X_1 = x_0 + x_1 * w   + x_2 * w^2
X_2 = x_0 + x_1 * w^2 + x_2 * w^4

Which is the same as:
X_0 = x_0 + x_1       + x_2
X_1 = x_0 + x_1 * w   + x_2 * w^2
X_2 = x_0 + x_1 * w^2 + x_2 * w

and it needs 4 (modular) multiplications.
However, apparently it can be also written as:
X_0 = x_0 + x_1           + x_2
X_1 = x_0 + x_1 * w       - x_2 - x_2 * w
X_2 = x_0 - x_1 - x_1 * w + x_2 * w

which uses only 2 multiplications.
I just cannot figure out how to derive - x_2 - x_2 * w out of x_2 * w^2, could someone please fill in the details for me? Are there other ways to reduce the number of multiplications? Is it possible to use just one (or even zero) multiplication? Does similar trick work for larger transforms?

What you have to show is that
$$W_3^2+W_3+1=0\tag{1}$$
with $W_3=e^{-j2\pi/3}$. It's a straightforward exercise to prove $(1)$.
In general we have
$$\sum_{k=0}^{N-1}W_N^k=0,\qquad N>1\tag{2}$$
with $W_N=e^{-j2\pi/N}$. Eq. $(2)$ can be shown by using the formula for the geometric sum.


Bayesian Information criterion is independent of prior - What does 'prior' mean here? I am trying to identify the possible number of states in my data. Each state corresponds to a different scenario. For example, I am measuring humidity and trying to identify the number of states in any given day. (Rain, Cloudy, Sunny).

For this I want to use the BIC to get the optimal number of states when modeling the data with a Gaussian Mixture Model. I am using scikit-learn for this and sklearn documentation says:

  Model selection concerns both the covariance type and the number of components in the model. In that case, AIC also provides the right result (not shown to save time), but BIC is better suited if the problem is to identify the right model. Unlike Bayesian procedures, such inferences are prior-free.

I am not sure what the priors are here.
Thanks in advance!

"Prior" refers to the "already known" probability (or its probability distribution function) of an event happening


What are the values on the x-axis when I plot the impulse response of a lowpass filter? I am trying to develop a filtering routine, using Java, and I am using the windowing method.

From my studies, I was told that if I make a plot of the impulse response of an ideal lowpass filter, I would get a SINC graph. In this situation, what am I plotting against (i.e., what values are on the X-Axis?)?

For digital filters the x-axis is samples.


Biquad coefficients using Magnitude (or Phase) Invariance Mapping Method I'm trying MIM (Magnitude Invariance Method) and PIM (Phase Invariance Method) for to improve biquad LPF response at low sampling rates. I'm looking some help and examples of usage if available.

EDIT_1:
I tried MIM method for 2nd order

LPF
fs=44100
fc=15000
Q=0.707

and

Peak filter
fs=44100
fc=10000
Q=0.707
db=6

Here are the resulting plots:

Responses looks promising.

Example code (LPF):

% Analog model for LP filter:
  fs = 44100;
  fc = 10000;
  Q  = 0.707;

  w0 = 2 * pi * fc;
  b2 = w0^2;
  a2 = [1 w0/Q w0^2];
  Ds = tf(a2, b2);

% Calculate coefficients using MIM
[Dz] = c2dn(Ds, fs, 'mim', 2, 4096*256);

bodemag(Ds, 'r', Dz, 'b');

2nd order MIM LPF results looks great for selected cutoff frequency (fs = 44.1kHz) as seen in plot (frequency range upto 23873Hz):

Coefficients for plotted 2nd order filter:

a [0.437417389844334 0.1220119196463959 -0.04071587087954469]
b [1 -0.6951293456643256 0.2138427842755107]

EDIT:

It looks like there are some problems in magnitude response when fs is increased and/or when fc is reduced. MIM/PIM paper mentions this as an issue.

So, it looks like MIM can be used for to improve 2nd order LPF magnitude response especially when using low fs and high enough fc. Magnitude issues can be bypassed by using higher order filters. Fortunately MIM method works correctly for my purposes.

EDIT: First and the second D(s) given in paper examples are wrong. The correct ones are (5s + 2) and (s + 10) / (s^2 + 4.1s + 0.4). With those all the MIM calculations become right. The third example (PIM) was correct.


What is the difference between the sampling frequency of signal and sampling frequency of filter I believe that there is no connection between the sampling frequency used for converting an analogue filter to digital filter and the one used to sample a signal that the filter will be used on. But I would like to confirm this from someone who knows.

I like to confirm that sampling frequency(e.g when using bilinear or impulse invariance method) used to convert analogue to digital filter can be chosen independent of the sampling frequency used to record analogue to digital signal. E.g. can an analogue filter converted to its digital equivalence using the bilinear transform at 100 Hz be used to filter a digital signal sampled at 1000 Hz?

I sadly don't have the reputation to comment on Fat32's answer, but let me try to answer directly instead: the frequency response of a digital filter is always relative to the processing sampling rate.
So if you design a filter with cutoff at 100Hz and 1kHz sampling rate, then you are really designing for a normalized cutoff frequency $f/f_s$ of 1/10 and whatever the processing rate (ie. the sampling rate of the processed signal), the cutoff will always be 1/10 times the processing rate. If you design a similar filter with cutoff at 200Hz and 2kHz sampling rate, then you will get the same filter.
Now, if you want to filter a signal sampled at 1kHz and you want the cutoff to be 100Hz, then you need to design a filter with a normalized cutoff frequency of 1/10. If you use the same filter on a signal sampled at 2kHz, then the cutoff will be at 200Hz. If you want to keep the cutoff at 100Hz for the 2kHz signal as well, then you need a different filter, this time for a (normalized) frequency of 1/20 instead.
The filter coefficients don't carry any information about the sampling rate, but rather everything is always relative to the sampling rate of the signal you are filtering and this applies to all digital filters, whatever the design method and whether FIR and IIR.


Unit impulse vs Kronecker delta vs Dirac delta? Apparently these 3 terms are similar, but what is exactly difference between them and what is exact value of each of them at time zero,1 or infinity?

Kronecker Delta used mainly in context of discrete time signals is defined as $$\delta_{ij} $$ is 1 only when $i=j$ or $i-j=0$, $0$ otherwise. So alternatively it is also written as $\delta[i-j]$. This value is exactly 1 and finite (refer below for the distinction with Dirac Delta)

The unit impulse is just the kronecker Delta with $j=0$ hence we only refer to unit impulse with one parameter $\delta_i$. Since $j=0$, this is alternatively written as $\delta[i ]$.Hence is 1 at $i=0$, $0$ otherwise.

The continuous time unit impulse can be thought of derivative of continuous time unit step (modified to have a slightly less steep than the sudden step at $t=0$). Refer to Oppenheim's book on signals and systems Chapter 1.

The Dirac Delta is used in context of continuous time signals to define impulses. It is defined as having an infinitely small width and hence a large magnitude, however the area under the curve integrates to 1. Refer to Oppenheim's book on signal and systems for various shapes that a Dirac Delta function could take, in all cases the area being 1.


fir filter - mean delay and snr from difference equation? I'm trying to get my dsp sea legs a bit, and am trying to complete a problem that asks for the mean delay and expected SNR boost for a given difference equation: y[n] = (x[n] + ... + x[n-N+1])/N

I'd love some general guidance on how to learn more about how to address this. Specifically,

I'd like to plot the filter function in octave/matlab, but am not sure how to get there from the difference equation. most of the examples I've come across assume that I have a frequency band to specify, but I'm not sure how I can derive that from just the difference equation.
I'm not sure how to calculate the expected SNR boost in the absence of any info about the noise source. Is it possible to calculate some concept of the amount of 'smoothing' the filter is expected to provide, and consider that to be the reduction in noise between the input and the output?

Thanks for any hints!

I'll give you some hints to help you solve the problem by yourself. Looking at the input-output relation, you can see that the output signal is computed only from a combination of delayed input signal values, and not from delayed output signal values. So your impulse response has a finite length and the filter can be implemented by a non-recursive structure (FIR filter). Using the impulse response, you can write the input-output relation as a finite convolution sum

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

Comparing (1) to your difference equation, you can easily see what the system's impulse response $h[n]$ is. Once you have the impulse response, you can compute the frequency response using an FFT. Plot its magnitude and phase to see how the filter behaves in the frequency domain.

The (group) delay is the negative derivative of the filter's phase response. You will find that the phase is a linear function, so the delay is a constant. You just need to determine the slope of the phase (either analytically or by inspection of the phase plot).

Without any further knowledge of your desired signal and the noise characteristics, you cannot say much about a possible SNR improvement.


How to implement a filter in matlab I want to implement the 'filter' function in matlab but I just can't seem to replicate the results I get when using the matlab function.

My understanding of the matlab function is that it takes 3 arguments: 'b' = array of numerator coefficients , 'a' = array of denominator coefficients and 'array' = the array to be filtered.So if you wanted to filter an array with a filter H(z) where:

$$
  H(z) = \frac{b_1 + b_2z^{-1}}{a_1 + a_2z^{-1}}
$$

and store it in an array called "result" then you would write

result = filter( [b1 b2] , [a1 a2] , array );

Now I tried to implement the 'filter' function ( with the number of coefficients in the numerator and denominator limited to 2) by getting the difference equation from H(z) then using it to get the resultant array but this route does not seem to work, ie

$$
  H(z) = \frac{Y(z)}{X(z)} = \frac{b_1 + b_2z^{-1}}{a_1 + a_2z^{-1}}
$$

$$
  X(z) (b_1 + b_2z^{-1}) = Y(z)(a_1 + a_2z^{-1})
$$

$$
b_1x[n] + b_2x[n-1] = a_1y[n] + a_2y[n-1]
$$

$$
\implies y[n] = \frac{b_1}{a_1}x[n] + \frac{b_2}{a_1}x[n-1] - \frac{a_2}{a_1}y[n-1]
$$

Meaning the matlab code should look something like

function [ result ] = myFilter( b , a , audio )
%MyFilter
%a,b = Two Element Arrays
%array = the array we want to filter

N = size(audio,2);

for i = 1:N
   if( i == 1 )
       result(i) = ( b(2)/a(1)) * audio(i);

   else
       result(i) = ( b(1) / a(1) ) * audio(i) + ...
                   ( b(2) / a(1) ) * audio(i-1) - ...
                   ( a(2) / a(1) ) * result(i-1);
   end
end

end

The function above does not produce the same results as the matlab filter function.Does anyone know what I am doing wrong? Have I understood how the filter function in matlab works?

You have a typo in:

result(i) = ( b(2)/a(1)) * audio(i);

That should be:

result(i) = ( b(1)/a(1)) * audio(i);


Angle modulation spectrum Is angle modulation spectrum always discrete or this property only holds for periodic modulating signal?

You will see lines in the spectrum of a phase-modulated signal if there is some periodicity in the modulating signal. If the modulating signal is not periodic, then the power spectrum of the modulated signal will be smooth.

As an example: consider a BPSK signal:

If the bits used to modulate the waveform are equally likely to be a one or a zero and are uncorrelated with one another, then there is no periodicity in the modulating signal. The power spectrum of a BPSK signal with independent equally-likely bits has a $sinc^2$ shape with no discrete lines.
If you then changed the data bits to add some repetitive framing data, for instance, causing a pattern of invariant bits in the modulating signal, you will begin to see lines appear due to the deterministic periodicity in the modulated signal.

In practice, one often makes an effort to ensure that bits are equally distributed to ones and zeros without repetitive patterns. Adding a source encoder (e.g. data compression) to your system can help to remove any repetitiveness in your data and whiten the transmitted bitstream.

One caveat to the above, courtesy of Dilip Sarwate's comment below, is duplicated here to make it more visible:

  However, modulating signals such as $\cos(\omega_1 t)+\cos(\omega_2 t)$ where $\omega_1 / \omega_2$ is an irrational number (i.e. no harmonic relationship exists) are not periodic but have lines in the spectrum at $\omega_c \pm m \omega_1$, $\omega_c \pm m \omega_2$, and $\omega_c \pm m \omega_1 \pm n \omega_2$ plus a lot of clutter (smooth spectrum) in between the lines. Here $\omega_c$ denotes the carrier frequency. There are also similar lines (and clutter) near the harmonics of the carrier frequency.


Are there paper about black box testing of optic flow algorithms? Do papers exist about comparisons of the resulting optic flow and its precision compared to the used optic flow method?

I would like to find paper that dealt with this topic in some way and had discussed different ways to conquer this problem.

Thank you in advance.

There is Middlebury, which evaluates almost all the algorithms and their evaluation results are online:

http://vision.middlebury.edu/flow/eval/results/results-n1.php

Or simply click any of the links at: http://vision.middlebury.edu/flow/eval/

They also have a paper on their research and evaluation methodology.


How to Save FFT output for fast read-in? For a project, I'm reading in hundreds of short (4 sec) wavefiles and giving them to scipy fft for the frequency information for further processing/experimenting in various ways. Problem is, I want to run the program many times while changing params, code, test new ideas, etc, but each time I'm rereading, re-fft-ing the original wavefiles. I'm wondering if there is a smart way to "save" the fft numpy arrays for quick retrieval, so I'd only have to read them in once, if this is possible and not a dumb idea.

That doesn't seem worth the bother.

FFT of a 7 second long wave file on my windows 10 laptop using Matlab takes about 5 milliseconds. So unless you have a particularly slow setup, you will not save a significant amount of time by storing the FFTs instead of just reading the raw wave file and doing the processing again each time.

That doesn't seem to justify the extra development effort for a non-standard file format.


How to identify impulse response of noisy system? My lecture slides mention that two-level ("2-tasoiset" in Finnish) signals are good in identification with linear systems. I am trying to understand this in noisy system. I cannot find any method really to identify the impulse response, perhaps good pre-processing and some other checks?

This is a sub-problem 2F here for the course Mat-2.4129 in Aalto University.

If a linear time-invariant system has (deterministic) input $x(t)$ and output $y(t)$, then the cross-correlation function $R_{y,x}$ of the output and input has value $h \star R_{x,x}$ where $R_{x,x}$ is the autocorrelation function of the input signal and $h$ is the impulse response of the linear system.  For completeness, the autocorrelation function $R_{y,y}$ of the output is given by
$$R_{y,y} = h\star \tilde{h}\star R_{x,x} = (h\star\tilde{h})\star R_{x,x} = R_{h,h}\star R_{x,x}$$
where $\tilde{h}(t) = h(-t)$ is the impulse response $h$ with time reversed, and $R_{h.h}$
is the autocorrelation function of the impulse response signal $h$.  In the frequency domain, this relation is the familiar power spectral density relationship
$$S_y(f) = |H(f)|^2 S_x(f).$$

Coming back to the business at hand, suppose that $x$ is a periodic signal with
a very long period, much longer than the measurable duration of the impulse
response $h$. Then, the output signal $y$ also has the same period, and essentially looks like periodic repetitions of $h\star \hat{x}$ where $\hat{x}$ denotes one period
of $x$..  The input autocorrelation function
$R_{x,x}$ is also a periodic signal and thus $R_{y,x}$ looks like periodic
repetitions of $h\star \hat{R}_{x,x}$ where $\hat{R}_{x,x}$ is the periodic
autocorrelation function of $\hat{x}$.

Finally, suppose that $\hat{R}_{x,x}$ is an impulse and $R_{x,x}$ is
thus a periodic impulse train with very long period.  Then,  $R_{y,x}$
is just periodic repetitions of $h$, the impulse response that we seek.
Now, this is all fine and dandy but everyone knows that impulses
don't exist in real life and money can't buy happiness. But as
the vulgar rich know, and all others suspect, money
can buy the most remarkable substitutes for happiness. Similarly, there exist
signals whose periodic autocorrelation functions approximate impulses.
If $\hat{R}_{x,x}(0)$ is very large while $\hat{R}_{x,x}(\tau)$
is very small, possibly even $0$, for all other $\tau$,
then $\hat{R}_{x,x}$ is essentially impulse-like. Even better,
such signals can even be chosen to be two-valued. They are, of course,
the binary $m$-sequences, also known as maximal-length linear feedback
shift register sequences. A $m$-bit linear feedback shift register
whose feedback connections are specified by a primitive binary
polynomial of degree $m$ produces a sequence of period $2^m-1$.
If the clock cycle duration is $T$ seconds, the output waveform
from the shift register has period $(2^m-1)T$ seconds. Regarded
as a two-level waveform with levels $\pm 1$, the periodic autocorrelation
function has a peak value of $(2^m-1)T$ which is very large
compared to the off-peak value of $-T$, and indeed the difference
can be made very large even with moderate values of $m$.

So, how does all this work? The input is a periodic two-level $m$-sequence
signal whose periodic autocorrelation is the proverbial inverted
thumbtack function. We compute the cross-correlation of the output
and input signals, effectively getting the impulse response of the
system as the result. True, the result is actually the response to
periodic thumbtacks but the discrepancy can be made small. Note also
that we get the impulse response without having an impulse
applied to the system. This is important: most linear systems
are actually nonlinear systems operated in a regime where the
signals are small in amplitude, and the excursions from the
nominal operating point are small enough that the tangent to
the (nonlinear) operating curve works well enough as a substitute
for the exact response (small-signal analysis). So, the input
is low-level but by a miracle of modern mathematics, we are
able to compute the response of the linearized version
of the nonlinear system to an impulse without using a large
signal that will blow fuses and drive the system into saturation,
overload, etc.


what is difference between Alias vs. down sampling? could someone help me explain downsampling in the context of Nyquist–Shannon sampling theorem and tell the difference between Aliasing and downsampling.Also is downsampling and undersampling the same ?Thanks in advance :)

"Alias" is to "downsampling" what "a broken leg" is to "jumping from a 40-feet height". An effect from a cause, with a very probable risk in general, but could be harmless if you are lucky.

I am not sure there is a clear shared difference between uses of downsampling and undersampling. To me, used as a verb, "to downsample" may suggest a rate reduction by an integer or a fraction. While "to undersample" could be more general (at some sample rate below its Nyquist rate). Some definitions equate undersampling and bandpass sampling.

Used as an adjective as in "the signal is downsampled" or "the signal is undersampled", I understand different meaning. In the first one, I suppose a voluntary action. In the second one, this can be interpreted as an observation: the signal has been sampled at a rate lower than expected.


Combining Analog Video Signals for Diversity Receiver System The system I want to build takes an analog video from two receivers (that are receiving from the same transmitter) and then combines the signals using equal gain combining (or maximal gain combining if a simple circuit exists to variably attenuator each signal).
My question is, how is this method practically implemented?
The best I've been able to find is that each analog signal is converted to digital and then has the gain adjusted before combining.

Is the combining of the signal usually done before or after the receiver?
Would the method of converting to digital and scaling gain introduce too much noise to the signal in practice?
What kind of ADC and variable amplifier would be needed in this application (NTSC composite video)?

Thanks for any help you can provide.

(or maximal gain combining if a simple circuit exists to variably attenuator each signal).

Nobody does maximum gain combining in analog circuitry (at least as far as I know). That would require something to estimate your channel coefficient, calculate its inverse, and convert it to a complex gain for an analog element - after having probably already converted the received signals to digital signals, to even estimate the channel coefficient.
You simply do that math in your digital signal processing. A single multiplication doesn't "hurt" much.
Equal gain combining is indeed sometimes done in analog hardware, but it's then usually called "phase shifting for analog beamforming" and mostly done in analog, because the realization in digital domain is too expensive due to a very high signal bandwidth or antenna count, typically in radar systems. You're not in one of these use cases, far as I can tell.

My question is, how is this method practically implemented?

Honest answer: very rarely... the era of analog video transmissions ended before MIMO / diversity combining became a thing.
That can be explained, and it also hints at what you need to do when implementing this!
For any combining method to be useful, you need your video signal to take different paths to different antennas (otherwise: no diversity gain).
Now, since most video transmission things back when analog video was the dominant method were television, i.e. were done in an outside environment, there was a slight problem: The transmitter would typically be a tower that can be received with a mostly straight path by antennas on rooftops.
That doesn't really lend itself to any temporal diversity: that channel typically doesn't change over time.
Now, if there's more than one independent paths for multiple antennas attached to the same system, that typically means there's some big reflector (typically: mountain behind your receiver), and these paths overlay.
That introduces multipath into every antenna, typically: Let's say your analog TV signal was 6 MHz wide: at $c_0=3\cdot 10^8\,\frac{\text{m}}{\text{s}}$, that basically means that as son as the delay between different important paths comes into to the region of 50 m, you'll have a visible problem. For larger delays, it becomes such an important problem that it got a name: Ghosting.

source
You can see the echo!
Now, your combining system relies on the existence of multipath, effectively, because NTSC runs at frequencies and over distances where the antennas get to large to have multiple independent antennas in an indoor / handheld device. But if multipath occurs, this is what it looks like... Before gaining anything from diversity combining, in a real-world system, you'd first solve the multipath issue, because that's a severe problem.
Solving that effectively requires an equalizer. So you build an equalizer, which requires your receiver to know something about the transmitted signal, and then use that knowledge to estimate what a filter would look like that "reverses" the multipath.
Proving that this can be done is an early 1990's thing. Realizing that doing that on analog transmissions is more hassle than moving to digital signalling happened exactly at that point. So, to my knowledge, not a single analog video system with diversity combining ever hit a market.

Is the combining of the signal usually done before or after the receiver?

But now you know what you'd need to do: first cancel echos, i.e. have an equalizer, for each antenna. Your combining techniques don't (or at least, not in their vanilla form) work at all if the channel is more than a complex coefficient – and a multipath channel is a complex impulse response, so you need to solve that first. After you've built an equalizer for an analog signal, you can then combine them.
Note that it's really a pain to equalize analog signals, especially things that aren't well-shaped. And NTSC really isn't a great candidate for that, there's very little that you can rely on to give you a channel estimate. That means your analog transmission needs to be slightly extended to include some kind of channel sounding sequence - and that essentially means you're adding something from a digital world (a known preamble) into your analog signal.
Let's assume you do that. (hint: TVs with such techniques did exist... you likely didn't have one. Early on, that was probably a "premium" feature, and when it became affordable, analog TV was on its way out.)

Would the method of converting to digital and scaling gain introduce too much noise to the signal in practice?

You do realize any modern transceiver, be it your digital TV, your smartphone, your wifi access point or your bluetooth thermometer, does that? The noise introduced by an 8 bit ADC will be laughably low compared to the noise you got in analog TV reception, anyways.
This isn't the 1980s...

What kind of ADC and variable amplifier would be needed in this application (NTSC composite video)?

Now, this is a very hard question. The answer to that is: figure out the bandwidth of your signal, and build a complete receiver frontend (which is a pretty complicated thing at first) that can digitize that bandwidth (Nyquist tells you how much sampling rate you need; if this is news to you, you'll have to read up on digital signal processing first). But that's just a lower limit: what ADC you pick depends on how much effort you want to spend on analog filters, what your system's digital interfacing looks like, what skills you have with FPGAs (you can't just plug an ADC into your computer)....
Effort-wise, and probably even financially, you'd just buy an SDR device instead. In these bandwidth regions, things like hackRFs, bladeRFs, even LimeSDRs become pretty affordable; I personally like (but also: am affiliated with) the Ettus product line of USRPs better. You'll need two channels that are synchronously sampled.
Now, after that you'll need to invent a viable method of adding some channel-sounding sequence to your NTSC signal (standard NTSC won't do!), then you need design an equalizer in your digital signal processing, then you can throw your combining at things.
Uff. At that point, really, moving to digital video would have saved you a lot of effort: First of all, digital terrestrial TV systems solve the multipath issues, typically through multicarrier (specifically: OFDM) methods. That directly gives you something you can combine.
Now, look at the quality of images that NTSC offers: you get interleaved video with around 525 lines... After applying a very well-intending smoothing filter to make things look less ugly, something like this:

Honestly, for the same video quality, you'd have to use far, far less bandwidth if you threw digital video encoding at the problem, you'd get a digital receiver, which is way easier to deal with if you're equalizing, and it's kind of not technology that's not been dead as a doornail for 20 years; I'd drop the NTSC requirement, and plug together something narrower in spectrum using the DVB blocks that come with GNU Radio (which I'm also affiliated with).
That gives you digital transmitter whose output you can either plug into a transmitting SDR (see above) or you can use in a simulation (which I strongly recommend before dealing with RF, because that opens up a completely new box of problems, like non-reproducibility of channel realizations).
You can use GNU Radio then to simulate your different channels, and then combine the resulting signals, to your heart's delight. GNU Radio's DVB toolbox brings both the transmitting as well as the receiving side for digital video broadcasting standards (and arbitrary modifications of these standards), and you can omit e.g. the channel coding and decoding for your experiments.


VAD speech databases Does anybody know any open voice/speech annotated database for Voice Activity Detection? Or any other corpus for speech evaluation purposes?

The ones I know are private: AURORA, TIMIT, SpeechDat-Car, Speecon, Catalan Speecon.

The other option is to create hand-labeled annotations for own recordings, or use a method like the Short-Term Energy Contour. But in the first case the amount of files to process would be insurmountable, and in the second, the energy based methods only work fine for free-noise recordings.

Researchers from the Johns Hopkins University have recently released a corpus of music, speech, and noise which, according to them, is suitable for training models for voice activity detection and music/speech discrimination.

See https://arxiv.org/pdf/1510.08484.pdf for details.


Why are ideal filters not particularly useful/ hard to do? Why is a LPF ideal filter:

a) Not useful?

b) Hard to make?

I was asked to make a Low pass filter (simulink block) that will take an input of a cutoff frequency, and for a signal with a frequency higher then that cutoff frequency, the amplitude would be zero & with zero roll-off. Hence ideal. Is this an impossible task?

You can only make a finite-order digital filter, since you don't have infinite CPU capacity, infinite RAM, infinite numerical precision and infinite time to wait.

With a finite-order filter, the Z-transform will be an analytic function over the complex domain. The frequency response is the Z-transform evaluated over the closed path $z = e ^{j\Omega}$ (i.e. the unit circle), where $\Omega$ is the frequency.

An analytical function evaluated over a path will result in an infinitely derivable and continuous frequency response (property of analytical functions, as long as you don't step over poles).

An ideal LPF is discontinuous at the cutoff frequency. So not possible to implement.

In practice, since the ideal LPF is a discontinuous frequency response, adding more coefficients to your filter will converge slowly to the desired frequency response (convergence is in mean-quadratic sense, which is slow). So you will end up with a filter with many coefficients, bad frequency response, lot of delay, and lots of roundoff errors.


Signal Timing Recovery, Loop Bandwidth and SNR As explained in the answer here, for the symbol timing recovery loop we need to keep loop bandwidth sufficiently below the modulation bandwidth of the signal so as to avoid the filtering of the signal itself.
But when we send a random sequence of bits the corresponding baseband signal will contain frequencies from $-\frac{f_s}{2}$ to $+\frac{f_s}{2}$, where $f_s$ is the symbol rate.
Since the signal has spectrum even at low frequencies, how can we ensure that the signal will not be filtered by the timing recovery loop?
I am asking because of the following reason:
1. I generate the PSK modulated data and add phase noise on it

dataMod = pskMod(dataVec, 8, 'gray');
tx_out = dataMod.*exp(1i*ph_out);

Here, the ph_out is the phase jitter signal and dataMod is the PSK modulated signal of the dataVec symbol stream.
2. I pass the phase of the data to the carrier recovery loop:

pll_out = pll(angle(tx_out), 0.005, 1e-6, 1e-3, 1);
pll_out1 = pll(ph_out, 0.005, 1e-6, 4e-3, 1);
pll_out2 = pll(angle(dataMod), 0.005, 1e-6, 1e-3, 1);

I create three different outputs, pll_out is the loop output with the complete signal, pll_out1 is the loop output with only the phase noise as input and pll_out2 is the output with only the signal (without noise) as the input.
I see the following at the output of the tracking loop:
1. The power spectral densities for input and output phase noise:

From this, I can see that the low frequency phase noise is getting filtered as expected.

The power spectral densities of input and output signal:

From this, it is seen that the signal power at low frequencies is filtered. But because of this filtering additional "noise" or "distortion" (whatever is the correct term) causes the following PSD at the output (yellow):

The frequency roll-off is explained very well by Dan here. But somehow now, the roll-off is not the problem for me.
3. The total noise at the PLL loop output is thus:

So, although the signal distortion/noise does not affect me at high frequencies (and its roll-off doesn't matter anymore) but my low frequency noise gets degraded.
So, how can we ensure that most of the signal power is outside the PLL bandwidth?

Ultimately we are trading how much signal we lose versus how much noise (phase noise dominantly) we are rejecting along with being able to track dynamics in the overall model between trasnmitter and receiver. If I was optimizing this I would do a simulation with the actual parameters for all noise sources and expected dynamics in my system, but I start with the rule of thumb as I suggested in this other post Loop bandwidth for symbol timing recovery of using something between $R/100$ to $R/20$ where $R$ is the symbol rate. The bottom graph shows the effect of a lower cut-off frequency in that more of the symbol energy is maintained.

To see this consider the following graph showing a bi-phase signal (+1/-1) after passing through a high pass function. This also shows the motivation to not have an excessively long stream of all 1's or all 0's in the modulated data, which is one reason for data scrambling at the source (another reason is timing recovery).

The top graph shows the bi-phase signal with infinite bandwidth including DC. The issue with maintaining bandwidth to DC in the receiver is phase noise and timing drift. Ultimately the DC value that here is centered on 0 would drift as the receiver estimate of the transmit carrier shifts from true value. We use the received signal to track this, which is a high pass function. If we track too quickly, as in the second graph, we remove too much of the symbol energy since the carrier tracking loop with wider bandwidth would assume the new DC value representing the symbol is actually the carrier, thus removing it. A very long stream of 1's or 0's in a row would have the same effect.

So the answer is to understand first what is the longest possible run of symbols that would not change in the receiver and then choose a loop bandwidth that minimizes the amount of symbol energy lost over that time interval, while ensuring the bandwidth is fast enough to track changing conditions and minimizing noise from low frequency sources such as phase noise.

As a rough measure of the impact of a loop BW using a loop BW of $R/20$, where $R$ is the symbol rate, consider that the single sided bandwidth of a properly pulse shaped waveform is slightly over $R/2$. If the loop BW was $R/20$ the SNR loss would approximately be 0.46 dB from tracking loss alone:

$$10log_{10}\bigg(1-\frac{R/20}{R/2}\bigg) =10log_{10}(0.9) = -0.46 dB $$

This can be decreased by lowering the loop bandwidth further (-0.18 dB for $R/50$), but you need to then see if you are adding more noise to this from the increased phase noise that would be added from the local oscillator and other jitter sources (sampling clock) as I show here PLL for Phase Demodulation and Carrier Tracking
, or if the dynamics in the system (Doppler for a moving transmitter or receiver for example) are changing faster than this can track.

This graph I have shows the basic concept and to consider the optimization problem, and as depicted there typically is a long shallow bottom in choosing the carrier recovery loop bandwidth where the noise is minimized. In a properly designed receiver, amplified thermal noise would be the dominant noise source (the goal of analyzing the cascaded noise figure), and the carrier loop bandwidth should be within this region where that applies. LO and jitter noise sources are typically set to be 10 dB or more below this noise level so as to add only a small contribution to the overall noise figure budget. As we see in this graph, if we however set the loop bandwidth too low, the LO/Jitter contributed noise sources will start to dominate. If we set the loop bandwidth too high, then tracking noise (which is tracking out the modulation as I show above as well as self-noise that can be originate from our carrier tracking loop) will start to dominate. Overall this is a system design challenge as it considers all noise sources in both the analog and digital domain and is a great example of how such mixed signal implementations cannot be designed individually in a vacuum.


Can filter "depth" be adjusted by mixing dry and wet signals? Can filter "depth" be adjusted by mixing dry and wet signals?

I.e. can I simulate e.g. a +6dB bandshelf/peak filter at 1kHz by mixing in some of the dry unequalized signal and some of a wet signal that has been bandpass filtered at 1kHz and the filter has around the same shape as the bandshelf/peak.

Can it be theoretically the same?

adding the input to the output of a scaled 2nd-order bandpass IIR will get you the classic peak/cut EQ curve.

adding the input to the output of a scaled 1st-order LPF or HPF will get you a 1st-order low-shelf or high-shelf EQ.

adding the input to the output of a scaled 2nd-order LPF or HPF will get you a particular form of a 2nd-order low or high shelving EQ, but there are many others.  doing it this way will often leave an unintended null or lump in the frequency response.  if you want a symmetric (in log frequency) and perfectly monotonic gain vs. frequency curve, i might suggest referring to the Audio EQ Cookbook.


RMS Amplitudes of voltage I am looking here at the course of current fluctuations and have calculated the RMS values, but somehow it seems strange to me: with the black line there are small fluctuations, but the RMS value is still comparably large with, for example, the blue...
The RMS values are:
red:0.83mV
blue:0.44
green:0.52
black: 0.31mV
compare with the picute

RMS includes the mean. If you want a metric for "variation" use the standard deviation instead or subtract the mean first (which is basically the same thing).
For the black curves moves slowly, but it still moves by a good amount. If you want a metric of the "fuzziness" than you can apply a high pass filter first that gets rid of the slow "drift" movemenent.


Fixing integer overflow in wave file I have this recording from my DAW that showed a strange form of distortion: the signal clipped at many places, but instead of simply being clipped (loss of information), somehow the integer value of the signal overflowed and has been stored in the file. It resulted in worse-sounding distortion, but the information is still there, as a mirrored image in the opposite side of the signal.
This is plainly evident when looking at the waveform. You can clearly tell when the signal clips because the rest of the wave is visible between the peaks.
How should I approach fixing this file? I am very comfortable with SOX, ffmpeg, the C language, SuperCollider and various other tools.

What you have there is "wrap around" and a bad case at that.
It's easy enough to fix though. Let's call your samples $x[n]$ and your max amplitude $x_{max}$, so your signal is currently bounded to $[-x_{max},+x_{max}]$
Go through the waveform sample by sample. Whenever you see a jump that's larger than your max amplitude, add or subtract twice the max amplitude to bring it back to where it was.
Something like
for (i = 1; i<n; i++)
{
   sampleDiff = x[i]-x[i-1];
   if (sampleDiff > xmax) { x[i] = x[i] - 2*xmax; }
   if (sampleDiff < -xmax) { x[i] = x[i] + 2*xmax;}
}

Make sure you are using a data type that accommodate a data range of at least  $[-2x_{max},+2x_{max}]$, which is quite likely NOT the original data type. Also make sure you don't start in a "wrapped" region, i.e. x[0] needs to be good


Simple Impulse Response I am trying to find the impulse response to this system:

At first it looks easy enough, $y(n)=x(n)(-1)^n$, set $x(n)=\delta(n)$, which makes $h(n)=(-1)^n$. But something seems wrong.

What am I missing?

(considering that $(-1)^n=\cos(\pi n)=e^{j\pi n}$, this system seems more complicated than it looks. It seems like $y(n)$ is $x(n)$ with a 180$^\circ$ phase shift.)

as the comments to your question mention you have not done the first step: system classification

Is the system linear? We have to input three different signal in the system and compare the outputs: the two signals $x_1$ and $x_2$ and the sum signal $x_{1+2} = x_1 + x_2$

\begin{align}
 y_1 &= (-1)^t x_1 \\
 y_2 &= (-1)^t x_2 \\
 y_{1+2} &= (-1)^t (x_1 + x_2) \\
\end{align}

and even at this point it is easy to see that $y_1 + y_2 = y_{1+2}$ that means the system is linear

Is the system time variant? it's usually very easy to say yes here when your system is a function of the time (or n in your case) $z(t) = (-1)^t $ Mathematically you have to input a single signal $x(t)$ at two points in time and check if the system gives you different outputs at different times:

\begin{align}
 y_1 &= (-1)^t x(t) \\
 y_2 &= (-1)^t x(t-\Delta t) \\
\end{align}

to compare 1 and 2 you have to time shift the output of the first system to make sure they overlap the way you want.

\begin{align}
 y_1 &= (-1)^t x(t) \\
 y_{1-\Delta t} &= (-1)^{t - \Delta t} x(t - \Delta t) \\
 y_2 &= (-1)^t x(t-\Delta t) \\
\end{align}

the comparison shows that those two are not the same, which means your system changes over time

\begin{align}
(-1)^{t - \Delta t} x(t - \Delta t) &= (-1)^t x(t-\Delta t) \\
(-1)^{t - \Delta t} &= (-1)^t  \\
(-1)^{ - \Delta t} &= 1  \\
\end{align}

this is not true in general, that means your system is linear and time variant and you can not use the description of a LTI System. But there's a definition of the impulse response of a time variant system that you could use:

Time-varying "impulse response"

\begin{align}
y(t) = \int h(\tau, t) x(t - \tau) d\tau
\end{align}

since your system is memory less and kausal (i don't think there's a trivial mathematical proof for those two properties but it is easy to see), your time variant impulse response is rather simple:

\begin{align}
h(\tau, t) = \delta(\tau) (-1)^t
\end{align}

which leads to the collapse of the integral in the time variant convolution (just check the rules for integration over the dirac function if you don't remember them):

\begin{align}
y(t) = \int \delta(\tau) (-1)^t x(t - \tau) d\tau \\
y(t) = (-1)^t x(t)
\end{align}

And the last line ( and basically your first line) is the best description for linear memory less time variant systems - a simple multiplication with a time variant function or in general just:

\begin{align}
y(t) = x(t) z(t) \\
\end{align}

bye bye stackexchange you managed to annoy me enough in my first two answers for me to decide not to come back.


DFT zero-padding of signals starting before n=0 If a signal starts before n=0, what part of the signal should be used to compute DFT after zero-padding? For example, x(n) = {1, 2, 3, 4, 5}, where x(-2) = 1 and x(0) = 3. If this signal is zero-padded to N=8, the new signal would be x'(n) = {1, 2, 3, 4, 5, 0, 0, 0}, where, again, x'(-2) = 1 and x'(0)=3. However, to compute the DFT of this zero-padded signal x'(0), should I use {1, 2, 3, 4, 5, 0, 0 , 0}, which starts with x'(-2), or should I use {3, 4, 5, 0, 0, 0, 1, 2}, which starts with x'(0) with the first few terms (n<0) folded over? In terms of the DFT X(k), I think its magnitude would be the same for the two approaches, but its phase would be different, right?

Use the second one Tony... It yields the correct implied phase relationship.


How to get RMS to Frequency chart for discrete samples of acceleration? I’m reading a book on Motorcycle dynamics and want to compare the vibration profile of my motorbike against this chart:

I am not sure how to get the RMS of acceleration at different frequencies.

I have a mobile phone with an accelerometer that samples acceleration at 100 Hz. I get filtered values of vertical acceleration from external sources (no gravity component).

I can do a Fast Fourier Transform (FFT) of 512 samples of data to get 256 amplitudes of waveforms for a short window of collected data. FFT contains frequencies 1-50Hz

Or I can square each accelerometer output, making them all positive, add them together and take a square root to get a single RMS value. In this case I’m not sure how to get the frequency component of these values.

Can someone help me understand how I get the chart like above from a continuous set of accelerometer output?

Here’s my attempt at plotting the log of amplitude of the FFT (in blue 0-50Hz). I don’t see any dominant component within the FFT output: at the same time the yellow chart above registers acceleration spikes over 0.155(g) threshold. I can see individual impacts on the suspension within the yellow chart, but am not sure how to get their frequency (impacts per second)

I also see this suggestion, but the octave example is not very clear. Is this relevant to my question?

If your time domain data is in units of vertical acceleration, then you can simply multiply the normalized (divide by N) raw fft values times the complex conjugate of itself and then take the square root. Each result will be the rms acceleration at that frequency.

Note, if your data is sampled at 100 Hz, you are going to see the content from 0 to 50 Hz of all frequencies that are within the analog bandwidth of your sampler. Any higher analog frequencies will "alias" into the 0 to 50 Hz band. I am not sure what you have for an anti-alias filter, but it certainly can't be a "brick-wall" filter right at 50 Hz. A properly designed filter may have hope of sufficiently rejecting signals above 60 Hz, in which case you can trust your data from 0 to 40 Hz (for example; I am more illustrating the need to understand this before interpreting your data).

Regarding "Normalized FFT"; some FFT implementations will have this built in while others you will need to normalize. You can test this by taking an FFT of 512 ones [1 1 1 1 ...1] and see if the first sample of the FFT is 512 or 1. If 512 then you need to divide the FFT result by the number of samples (512) to normalize the result.  And by "raw FFT" data I am referring to the 512 complex FFT samples, not the FFT frequency spectrum you plotted above which I assume is the power spectrum which already performed the complex conjugate multiplication and is only showing the positive frequency axis, likely with a dB vertical scale.

Usually these limits versus frequency are also normalized to be per Hz bandwidth. The FFT noise bandwidth is 1 bin (without further windowing which is another topic of importance if you start seeing dominant tones). So if you sample at 100 Hz with 512 bins, then each FFT bin has a noise bandwidth of $100/512= 0.1953$ Hz.  A more accurate measurement would average the power (the complex conjugate multiplication of the normalized FFT) over 5.12 bins before taking the square-root to get the average rms value per Hz bandwidth. A simpler and more conservative approach (will bound the actual max limit) would be to multiply the power in each bin by 5.12 before taking the square-root.

You may see spikes above the threshold in the time domain, but the limit is rms and over frequency, so those results will expect to be less.


code construction process in polar coding Arikan in his paper (arXiv link or IEEE link) about polar coding, explains how you can force the channel to polarize by applying a linear transformation (the generator matrix). He calls this a split and combine process, where he creates synthetic channels that have near perfect or near useless capacity.
But I understand that this is happening so that you can find the best positions in your block to place information bits to be coded. So what does that have to do with the actual channel?

Polar codes exploit the channel polarization phenomenon. But it helps if you separate the channel polarization and polar coding for a while.
In channel polarization, the bits that you need to transmit i.e., $\mathbf{u}=[u_1,u_2,\dots,u_N]$, are transformed into $\mathbf{x}=[x_1,x_2,\dots,x_N]$. Each of the transformed bit $x_i$ is transmitted over the actual channel $W(y|x)$. From viewpoint of the transformed or coded bits, nothing special happens as each coded bit $x_i$ goes through the same channel and is recieved as $y_i$. However, looking from the view point of all the user or uncoded bits together, $\mathbf{u}$ goes through a vector or block channel and is received as $\mathbf{y}=[y_1,y_2,\dots,y_N]$. This is what is referred to channel combining. The combined channel, $W_N$, is then split into $N$ virtual/synthetic/artificial channels, one for each of the input bits $u_i$. Now this is where magic happens i.e., they way these virtual channels are defined. Virtual channel of each bit $u_i$ has the following outputs:

Output of the combined channel $\mathbf{y}$. Nothing special about this as this is what we get from the actual channel $W(y|x)$ as the received codeword.
Exact values of $u_1,\dots, u_{i-1}$. This is the special part. This implicitly assumes successive cancellation decoding. If the true knowledge of  $u_1,\dots, u_{i-1}$ is available, the virtual channels $W_N^i(\mathbf{y},u_1,\dots, u_{i-1}|u_i)$ exhibit polarized behavior.

Notice that nothing happened to the actual channel $w(y|x)$ through which the coded bits travel. It is like performing MRI with the use of contrast dye. You can do the imaging without the contrast dye. But if you inject the dye into your blood (channel combining) and then do the imaging (channel splitting), you get better results.
Now we come to polar coding. In code construction stage, we use channel polarization to figure out which bits to freeze. Then when the code is used for error correction, the encoder does the channel combining and the use of successive cancellation decoder actually does the channel splitting.


Understanding the FFT phase spectrum with a simple example I'm trying to compute the DFT using scipy's functions. I don't understand why the phase spectrum of a simple sine wave with 2 Hz frequency doesn't show $\pm\pi/2$ at the $\pm 2Hz$ frequencies. Instead, the phase plot seems to have some linear dependence in the frequency, which I don't understand. I provide the code for assistance. How can this be fixed? Looks like a simple issue I'm not grasping. Please help.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftfreq

# 1Hz sine wave
npts  = 100
tmax  = 10
t     = np.linspace(0, tmax, npts)
y     = np.sin(2*np.pi*2*t)
dt    = tmax/npts

# FFT computation
Y      = fft(y)
freq   = fftfreq(npts, d=dt)
amplit = abs(Y)/npts
phase  = np.angle(Y)
phase  = phase / np.pi

fig, ax = plt.subplots(1, 2, figsize=(10,4))
ax[0].plot(freq, amplit)
ax[1].plot(freq, phase)

ax[0].plot([2,2], [0,0.5], '--')
ax[0].plot([-2,-2], [0,0.5], '--')
ax[1].plot([2,2], [-1,1], '--')
ax[1].plot([-2,-2], [-1,1], '--')

For sinusoids that are not exactly integer periodic in the FFT length, an FFT measures the phase at a circular discontinuity.  And that discontinuity flips direction as frequency changes from slightly below to slightly above an exact integer periodic-in-aperture frequency.  This is part of the effect of the default rectangular windowing of any finite length DFT or FFT.  Thus making the FFT phase result of a frequency sweep hard to interpret.

The way to remove this measurement at a circular discontinuity is to shift the waveform discontinuity away from the phase measurement point (FFT input sample 0) by doing a (circular) FFTShift before the FFT.  Then your FFT phase result will be referenced to the center of the original input vector (original waveform sample N/2), which will be continuous for continuous input functions or waveforms.

After the FFT following an fftshift, you can use the FFT frequency and phase results for the the center of the original input waveform to compute the phase at any other point (for instance the beginning, or sample 0 of the original input waveform.  Or the end, etc.).


Gain of FIR/IIR filters What are the formulas for signal and noise power gain of digital filters (FIR and IIR)?

For a FIR, I've seen in Harris' windowing paper that the DC gain is the sum of the filter weights.

$$
G = \sum_{i=0}^{N-1} w_i
$$

For a FIR, I've seen that the noise gain is the square root of the sum of the square of the weights.
$$
G = (\sum_{k=0}^{N-1} w_i^2 )^\frac{1}{2}
$$

What about an IIR?

Is there a reason you don't think that the general formulae are:

$$
G = \sum_{k=0}^{\infty} h[k]
$$

and

$$
G = (\sum_{k=0}^{\infty} h^2[k] )^\frac{1}{2}
$$

These will work, provided the system is linear and time-invariant (and BIBO stable).

To work with a rational transfer function:
$$
H(z) = \frac{\displaystyle\sum_{m=0}^{M} b_m z^{-m}}{\displaystyle\sum_{n=0}^{N} a_n z^{-n}}
$$
then the DC gain (at $z = e^{j0}$) will be:
$$
G_{DC} = \left. H(z) \right|_{z=1 = e^{j0}} = \frac{\displaystyle\sum_{m=0}^{M} b_m}{\displaystyle\sum_{n=0}^{N} a_n }
$$

I'll have to think a bit more about the second gain in your question.


Calculation of gain at certain frequency in case of nth order IIR filter There is this method for to set 0dB gain to be at wanted frequency (fc) (Octave/Matlab example for biquad LPF):
% needed for Octave -------------------------
pkg load signal
% -------------------------------------------
clf;
% calculate coefficients --------------------
fs = 44100; % sample rate
fc = 700; %Hz
fpi = pi*fc;
wc = 2*fpi;
wc2 = wc*wc;
wc22 = 2*wc2;
k = wc/tan(fpi/fs);
k2 = k*k;
k22 = 2*k2;
wck2 = 2*wc*k;
tmpk = (k2+wc2+wck2);

a0 = 1;
a1 = (-k22+wc22)/tmpk;
a2 = (-wck2+k2+wc2)/tmpk;

b0 = (wc2)/tmpk;
b1 = (wc22)/tmpk;
b2 = (wc2)/tmpk;

b = [b0 b1 b2];
a = [a0 a1 a2];

FLT1 = tf(b, a, 1/fs);

% adjust 0dB @ 1kHz -----------------------------
fc = 1000;  % Hz
w = 2.0*pi*(fc/fs);
num = b0*b0+b1*b1+b2*b2+2.0*(b0*b1+b1*b2)*cos(w)+2.0*b0*b2*cos(2.0*w);
den = 1.0+a1*a1+a2*a2+2.0*(a1+a1*a2)*cos(w)+2.0*a2*cos(2.0*w);
G = sqrt(num/den);

b0 = b0/G;
b1 = b1/G;
b2 = b2/G;

b = [b0 b1 b2]
% ------------------------------------------------

FLT2 = tf(b, a, 1/fs);

% plot
nf = logspace(0, 5, fs/2);
figure(1);
[mag0, pha0] = bode(FLT1,2*pi*nf);
semilogx(nf, 20*log10(abs(mag0)), 'color', 'g', 'linewidth', 2, 'linestyle', '-');
hold on;
[mag, pha] = bode(FLT2,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'm', 'linewidth', 2, 'linestyle', '-');
legend('LPF', 'LPF 0dB@1kHz', 'location', 'southwest');
xlabel('Hz');ylabel('dB');
axis([1 fs/2 -30 15]);
grid on;

How are formulas for to resolve num and den derived so calculation of G for nth order filter can be done? As for an example for a 4th order filter:
a = [1.00000  -0.61847  -1.09281   0.43519   0.30006];
b = [6.9411e-03   1.1097e-02   5.2508e-03   6.9077e-04  -3.2936e-06];
fc = 1000;  % Hz
w = 2.0*pi*(fc/fs);
num = ...; % ????
den = ...; % ????
G = sqrt(num/den);
b(1) = b(1)/G;
b(2) = b(2)/G;
b(3) = b(3)/G;
b(4) = b(4)/G;
b(5) = b(5)/G;

You just need to evaluate the transfer function on the unit circle at the frequency of interest:
$$H(e^{j\omega_0})=\frac{\displaystyle\sum_{k=0}^Nb_ke^{-jk\omega_0}}{\displaystyle\sum_{k=0}^Na_ke^{-jk\omega_0}}\tag{1}$$
and take the magnitude.
For the special values $\omega_0=0$ and $\omega_0=\pi$, Eq. $(1)$ simplifies to
$$H(e^{j0})=\frac{\displaystyle\sum_{k=0}^Nb_k}{\displaystyle\sum_{k=0}^Na_k}\tag{2}$$
and
$$H(e^{j\pi})=\frac{\displaystyle\sum_{k=0}^N(-1)^kb_k}{\displaystyle\sum_{k=0}^N(-1)^ka_k}\tag{3}$$
respectively.
EDIT: If you want a formula that directly expresses the squared magnitude of $H(e^{j\omega})$ then use this:
$$\big|H(e^{j\omega})\big|^2=\frac{\displaystyle r_b[0]+2\sum_{k=1}^Nr_b[k]\cos(k\omega)}{\displaystyle r_a[0]+2\sum_{k=1}^Nr_a[k]\cos(k\omega)}\tag{4}$$
where $r_a[k]$ and $r_b[k]$ are the autocorrelations of the denominator and numerator coefficients, respectively:
$$r_a[k]=a[k]\star a[-k]\\r_b[k]=b[k]\star b[-k]$$
where $\star$ denotes convolution.


Embedded Correlation Function I'm trying to write a correlation function on an embedded device. Because its embedded, I need to make sure that the memory/CPU usage is minimum.

Here's the deal:
- I have a saved waveform sampled at 1 MHz
- I have timing info for the waveform
- I am going to be taking a new waveform at 200 kHz with new timing info
- The timing info for both waveforms is comparable, neither of them start at zero
- I want to determine the delay between one waveform and the next precise to 1 us.

The method I was going to use to do this was:
- Find the maximum of each waveform (there should only be 1 peak in each, and they should match up)
- Set the indexes of each array so that the maximums of the waveforms match up. I'll save the amount of clock cycles that I am shifting in order to do this.
- Then I will correlate this and shift a few samples in each direction and correlate as well. The highest correlation value corresponds to the correct amount of delay, and I will keep this value.

The problem is that my array manipulation skills are limited. I'm struggling when trying to implement this algorithm, especially because one of the waveforms is sampled at 5X the other. Does anybody have any references for this or knows how I can achieve this easily?

The specific problems I'm running into are:
- How do I set the indexes differently so that the maximums match up? I don't want to write a new register because this takes up memory and takes a long time in terms of read/write cycles.
- How do I compensate for the different sampling frequencies? I could fill the 200 kHz waveform with zeros (which won't affect the correlation) but this will take up more memory with no real purpose.
- The timing info isn't precise, in that they don't necessarily line up perfectly. For example, the timing numbers "10000, 10020, 10040, 10060, 10080, and 10100" for the 1MHz signal  doesn't necessarily match up with "10000, 10100" for the 200 kHz signal. So I can't rely on using the timing numbers as an index very easily.

If what I am asking isn't very clear, I can try to embellish a little on requests, but I don't want this to be too long.

How fast is the signal you are trying to analyze? The maximum theoretical frequency you can handle is going to be 100KHz (sampling rate / 2). Giving a "true" answer with 1 usec resolution is not possible with a 1MHz sampling rate. Practically this may not matter.

Without specific timing constraints it is hard to give a good answer. To get your signals to have the same sampling rate you can do linear interpolation between the points of your slow data or average your fast data.

I recommend averaging your fast data. If you are working in counts from an ADC summing 5 of the fast samples and dividing by 5 will be fast. Your compiler will likely use reciprocal multiplication to avoid the division operation.

Trying to line up the signals based on a single peak is likely to give incorrect results.
An alternative method is to maintain an auxiliary array with a selection of the N highest peaks stored. Do not simply get the N largest points; they will likely be clustered together and will not give a good result. Instead track the largest point and do not allow a new peak unless the signal has dropped below some threshold value (say 90% of the peak) to introduce hysteresis. The exact number of peaks should be chosen on based on the shape of your data and the desired level of confidence as well as memory and CPU time.

To match the peaks scan through the second data set maintaining a best fit offset. When you find a possible match to the first peak (a value within some tolerance) you can then compare the "future" peaks. Finding the offset with minimal error for the N peak points is not too onerous (sqrt((x_1 - xbar)^2 + ... (x_n - xbar)^2)).

This can be done in only one pass through each of the data sets with only the extra memory required being the auxiliary array that keeps the N peak points.


How to calculate the frequency response of LTI system based on measured impedance I have an analog circuit containing unknown resistance, capacitance, and inductance. I have measured its complex impedance across the frequency range of interest. Is there an equation that will allow me to compute the frequency response (complex voltage gain as a function of frequency) of the circuit given the impedance? Clearly lower impedance magnitude would result in greater voltage gain magnitude, but what is the exact mathematical relationship?

As an aside, I am aware that I could directly measure the frequency response e.g. by taking the Fourier transform of its impulse response. This question is not about how to measure the frequency response, it is about the mathematical relationship between the impedance and frequency response.

this is more about electronics, but i think it can be asked here.

when you are looking for a frequency response of an impedance, i believe this means you have a two-terminal (a.k.a. one-port) device.  so, instead of an input voltage and output voltage to relate by a transfer function (and from the transfer function you get a frequency response), what you have is a voltage applied to the two terminals and a current measured, or a current shoved into (and out of) the terminals and the resulting voltage measured.

the complex impedance of the two-terminal device, as a function of frequency, is the frequency response the device assuming current is the input and voltage is the output.  you can plot it on a dB scale, but without a reference resistance, you don't have a reference for 0 dB on that frequency response plot.  but you have relative dB for the different impedances at different frequencies.

maybe ask the EE guys this question.  i dunno.