Convolving time reversed discrete input with time reversed discrete impulse response

Here is another example of how incredibly poor choice of notation misleads students everywhere. I do not possess a copy of any of the various tomes on signal processing that seem to be revered as the fifth Gospel on this site, and so
I do not know if the book cited by the OP actually has written $(1)$,
but I assert that it is against mathematical common sense to write things like
$$y[n] = x[n] \star h[n].\tag{1}$$
It is generally accepted mathematical convention that is a symbol
has the same meaning everywhere it appears in an equation or expression.
For example, when we write
$$y(t) = x(t)h(t),~ -\infty < t < \infty, \tag{2}$$
we mean by this that for every choice of real number $t$, the value
of $y(t)$ is the same as the product of the values of $x(t)$ and $h(t)$;
for example, $y(3) = x(3)h(3)$, and $y(313.012) = x(313.012)h(313.012)$,
ans so on and so forth.
So how are we supposed to make sense of this monstrosity $(1)$?
Clearly $y[3]$ is not $x[3]\star h[3]$ whatever meaning we might
ascribe to the latter quantity; the value of $y[3]$ depends on
a lot of other values of $x$ and $h$.  The poor notation misleads
students into writing things like

$$\sum_{k=-\infty}^{\infty}{\left((-1)^{k}(u[-k]-u[-k-8])\right)\left(u[-n-k]-u[(-n-8)-k]\right)}$$

which is not the correct expression for the convolution (or the correlation,
for that matter) of
$$x[n]=(-1)^n{u[-n]-u[-n-8]}\quad \text{and} \quad
h[n]=u[n]-u[n-8]\tag{3}$$

The correct way to do this is write
$$y[n] = (x\star h)[n] = \sum_{k=-\infty}^\infty x[k]h[n-k]$$
and then substitute $n-k$ for $n$ in $h[n]=u[n]-u[n-8]$
to
get $h[n-k] = u[n-k]- u[n-k-8]$ and not $u[-n-k]-u[(-n-8)-k]$
the way that the OP has it.

Some initial thought might have
shown that (assuming that $u[n]$ is the unit step),
that $h[n]$ is $0$ for $n < 0$ (since then both $u[n]$ and $u[n-8]$
are $0$, and that $h[n] = 0$ for $n \geq 8$ (since then
both $u[n]$ and $u[n-8]$ are $1$. Consequently for any given
fixed integer $n$, $h[n-k]$ is nonzero only for
$$0 \leq n-k \leq 7 \Rightarrow n-7 \leq k \leq n$$
and so the limits on the sum could be reduced from
$-\infty$ and $\infty$ to $n-7$ and $n$ (cf. my comment
on the OP's question).  I will leave it to the OP to figure
out whether the fact that $x[n]$ is also time-limited
can be used to further reduce the limits on the sum.
Hint: the answer to this might depend on the specific
value of $n$, and you may need to try several different
values of $n$, e.g. $0, 3, 8, 10, -40$ etc to see if you
can discern a pattern.


Different definitions of unit step signal

The Wikipedia article on Heaviside has a section that discusses the 3 different conventions on $H(0)$:

https://en.wikipedia.org/wiki/Heaviside_step_function

Essentially, the step function is used in probability theory as well but is called something different.

I believe that the "symmetric" $H(0)=\frac{1}{2}$ version is probably most suitable in Linear Systems Theory. In discussions on the Gibbs phenomena, the truncated Fourier series interpolates discontinuities at the mid-value which in this case is $\frac{1}{2}$.

Also, if you look at the article, there are a number of limits of continuous functions like Erf, that approach the step as some parameter approaches some value. If you were to set up a physical experiment, you would never be able to realize an ideal step function, you would need infinite bandwidth, so in reality, the best you could achieve would be something like a one of those continuous functions.  $H(0)=1/2$ is a vestige of physical realizability.


Convolution in the time domain means you are multiplying their frequency transforms. This is one of the many properties of the Fourier Transform.

Now, when you flip a signal X[n], you can also write it down as x[n1], where n1 = -n. If you look at the time scaling property, try replacing 'a' by -1. You'll observe that the F(W) will now become F(-W).

Since F(-W) is exactly the same as F(W), your result does not change.


Something like that should do. Please keep in mind that this code can be further improved by windowing your signal, using decibel scale for spectrum, taking only first half of a spectrum, etc. This is all cosmetic stuff and you can do it on your own (or maybe I will edit this answer at some point).

%% Time domain parameters
fs = 20500;    % Sampling frequency
dt = 1/fs;     % Time resolution
T = 10;        % Signal duration
t = 0:dt:T-dt; % Total duration
N = length(t); % Number of time samples

%% Signal generation
f0_1 = 100; % fundamental frequency of first sinusoid
f0_2 = 200; % fundamental frequency of second sinusoid
x1 = sin(2*pi*f0_1*t); % first sinusoid
x2 = sin(2*pi*f0_2*t); % second sinusoid

% We want 200 Hz signal to last for half of a time, therefore zero-out
% second half - use the logical indexing of time variable
x2(t>5)=0;

% Now add two signals
x = x1+x2;

% Calculate spectrum
X = abs(fft(x))/N;
% Prepare frequency axis
freqs = (0:N-1).*(fs/N);

%% Plotting time and frequency domain
% Time domain signal plot
subplot(211)
plot(t, x)
grid on
xlabel('Time [s]')
ylabel('Amplitude')
title('Time domain signal')

% Spectrum plot
subplot(212)
plot(freqs, X)
grid on
xlim([0 max(freqs)])
xlabel('Frequency [Hz]')
ylabel('Amplitude')
title('Spectrum')

%% Playing back signal
% Normalize the audio:
x = 0.99*x/max(abs(x));

% For MATLAB R2014a use audioplayer
player = audioplayer(x, fs);
play(player)

% For older versions use wavplay
% wavplay(x, fs);

And apart from playing back it will also produce:


Estimating the PSD (Periodogram): How do you find the Expected Value

TLDR;

The mistake is here:

  If I calculate the expected value, which is according to
  definition an average, this way

  $E\left\{\widehat{P_x}(e^{jω})\right\} = \frac{1}{2π}\int_{<2π>}
> \widehat{P_x}(e^{jω})dω$

The mistake is in your interpretation of expectation. It is not an average over all frequencies.

For each fixed $\omega$, the expectation is taken over many realizations of the random signal $x$.

Gory details

Suppose you are given a complex valued random signal $x(n)$ for $0\leq n\leq N-1$. I would emphasize that this is a random signal. If you repeat your experiment you'll obtain a different set of $N$ values each time (called "realizations" of $x$).

Recall how the periodogram $\hat P_x(e^{j\omega})$ is computed. You first compute the sample autocorrelation at each lag $m$:
$$
\hat r_x(m) = \frac{1}{N-m}\sum_{n=0}^{N-m-1}x(n)x^*(n+m) \;\;\;\;\;\;\;\;\;\;\;\;(1)
$$
and then get the periodogram estimate by taking the DTFT:
$$
\hat P_x(e^{j\omega}) = \sum_{m=-(N-1)}^{N-1} \hat r_x(m) e^{j \omega m}.
$$
Since we started with a random $x$ we would like to know how this estimate performs on average i.e. average over the randomness in $x$:. So we take the expectation:
\begin{eqnarray}
\hat P_x(e^{j\omega}) &=& \mathbf E\left[\sum_{m=-(N-1)}^{N-1} \hat r_x(m) e^{j \omega m}\right]\\
&=& \sum_{m=-(N-1)}^{N-1} \mathbf E[ \hat r_x(m)] e^{j \omega m}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(2)\\
&=& \sum_{m=-(N-1)}^{N-1} \mathbf E\left[ \frac{1}{N-m}\sum_{n=0}^{N-m-1}x(n)x^*(n+m) \right ] e^{j \omega m} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(3)\\
&=& \sum_{m=-(N-1)}^{N-1} \frac{1}{N-m}\sum_{n=0}^{N-m-1}\underbrace{\mathbf E[x(n)x^*(n+m)]}_{\mbox{autocorrelation function of }x} \; e^{j \omega m}\\
&=& \sum_{m=-(N-1)}^{N-1} r_x(m) e^{j\omega m}  \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(4) \\
&=& P_x(e^{j\omega})    \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(5)
\end{eqnarray}

where (3) follows by substituting (1) in (2), (4) is by definition of autocorrelation $r_x$ of the random sequence $x$ and (5) is by definition of $P_x$.

Once again, note that the expectation was taken to "average out" the randomness introduced by $x$. It is not an average over frequencies because the frequencies are not random here.


Some doubts in the matlab code of computing Wigner distribution

The FFT comes from the summing operation and multiplication by the
complex pointer term in the formula. Compare it to the formula for
the DFT (or FFT):

$$X_k = \sum_{n=0}^{N-1} x_n\cdot e^{-\frac {2\pi j}{N}kn}$$
This is not similar, real just returns the real part of the array.
The Wigner distribution is real, but because of rounding errors, there is always a very small complex part in the result. This operation removes it and just returns the real part.


How can we measure the SINR of the nodes: as we do not have the Power x Gain product without the knowledge of Power and Gain of IoT nodes.

Not at all, you have zero information about the N in SINR, and without channel knowledge you also can't have any information on the S, nor the I, even if you knew the transmit powers (which you don't).

  Is it right to assume that the receiver of the AP won't be able to get any legible data (information) unless the SINR is above a certain threshold?

Yes. If you're into NOMA, you should be aware of a bit of basic information theory (otherwise, it's hard to understand why you'd even ever want to go non-orthogonal):

In this case, a very basic corner case of the Shannon Limit holds true. Some call it the ultimate Shannon Limit, but what it really just says is that if your bit energy is below -1.59 dB of the energy of your redundancy, you can't get any information from your reception.

In a multi-user system, one can (often iteratively) gain information about the interference¹, but if noise² alone pushes the SINR below -1.59 dB, you've lost and better redesign your system.

 ¹Thus reducing restoring the receiver's ability to "remove" parts of the entropy the interferer adds to the received signal.

²"Noise" actually being all entropy in the reception that doesn't come from the information source. That incorporates all entropy of interference that you can't "cancel", however that happens.


First, you should use informative variable names. For example, use "upsample_ratio" instead of x, "inSignal" instead of y etc. This will make your code much easier to read and understand.

Second, you should avoid using for loops in matlab if possible, and especially don't use nested for loops, since they take very long time in matlab.

Third, you overuse the variable y. Use intermediate variables to keep temporary and result values.

Forth, your algorithm concept is wrong in case the upsampling ratio is not a whole number.

Specifically for your code, start by replacing the inner loops with array operations:

vec = length(y):-1:i+1;
y(vec+q) = y(vec);
y(i+1:1:i+1) = 0;


1) Matching lines across multiple views is a common research problem and is reasonably well studied:
If you have the end points for example : http://cmp.felk.cvut.cz/~werner/software/lmatch/lmatch_memo.pdf

Even tough this deals with more geometric properties, if your scenes are well conditioned, you could as well use the gradient/intensity information to match the lines.

Here is a full source code on multiview line matching from monocular setups:
https://github.com/manhofer/Line3Dpp

You might also like to start with that.

2) Using a KALMAN filter on the lines would be a good idea. There are also more sophisticated algorithms. For example, in each frame you might want to minimize cost function to obtain the subpixel line, which has better resemblence to the line in the previous frame.

You can also track the camera pose and from there hypothesize the lines.

You try to avoid SLAM, but maybe, triangulating the lines and obtaining 3D structure would be the best in reducing the errors.


Multiplication and Convolution in Wireless Communications

The short answer is that in the first case the channel delay spread spans one symbol time duration or less (which is called frequency flat fading channel), and can be approximated by a complex number according to the central limit theorem, while in the second, the channel delay spread spans several symbol duration (which is called ISI or frequency selective channel).

EDIT1: Mathematically, the flat fading channel can be represented as $$h(t)=h_k,\,\,\,\text{ for } kT_s\leq t<(k+1)T_s$$ while in the frequency selective case it can be represented as $$h(t)=\sum_{n=0}^{N-1}h_n\,\delta(t-nT_s)$$ where $N$ is the number of symbols the delay spread of the channel spans, and $T_s$ is the symbol time duration.

Now suppose the transmitted signal is $s(t)$. The noise-free received signal for the first channel is $$r(t)=h(t)\star s(t)$$ which when sampled at time $kT_s$ for the $k$th symbol equals to $$r_k=h_k\,s_k$$ while in the second case it's $$r(t)=\sum_{n=0}^{N-1}h_n\,s(t-nT_s)$$ which when sampled at $t=kT_s$ equals to $$r_k=\sum_{n=0}^{N-1}h_n s_{k-n}$$ where for brevity I used the notation $s(kT_s)=s_k$.


What does it mean for a function to have frequencies?

I understand that "to have frequencies" can be misleading. This is a shorthand for "having non-zero energy or amplitude at two (specific) frequencies". Rephrased, for the cosine: when a  cosine signal is represented in the frequency domain, it has two non-zero components at frequencies $\pm \omega$, and the other values anywhere on the frequency axis have zero amplitude. Or: its frequency spectrum has only two non-zero values.

Frequency in signal processing often refers to a domain of analysis well-suited to certain signals and systems (linear systems for instance). With Fourier analysis, such signals can be represented as a weighted sum of sines or cosines, and essentially the frequencies a signal has is directly  related to the number of sines/cosines with non-zero weight in the sum. One could be more precise, depending on which type of signal and Fourier transformation one is talking about (discrete/continuous).

The concept of frequency is somehow independent on the signal, just like the time, and time and frequency are interlinked via Fourier transforms. So "having frequencies" by itself has no meaning for a signal, just as "having time". But a signal can be described with either amplitudes in the time domain, or amplitudes in a frequency domain.

In that sense, whenever a Fourier transform is properly defined for a signal (integrability, etc.), and there is at least one non-zero value, the signal "has frequencies".

So, most signals you will meet have frequencies, except the zero signal, and those whose Fourier transform cannot be defined correctly. As you might see later in signal processing lectures, one often uses "having frequencies" in a loose sense, meaning "above a certain threshold" instead of being simply non-zero, as in DSP practice, it is rare to have exactly zero values (finite signals, approximations, etc.)


I found that the reason for the spike on the right is due to the inclusion of pulse width of time along the time axis.


How to differentiate a time domain signal in the complex transfer function?

The transfer function is commonly defined as output divided by input, and not so often the other way around (as you did in your question). Also Matlab uses that more common definition. So the transfer function $H(s)$ between current $I(s)$ (input) and velocity $V(s)$ (output) would be defined by

$$H(s)=\frac{V(s)}{I(s)}\tag{1}$$

Since acceleration $A(s)$ is the derivative of velocity, we have

$$A(s)=sV(s)\tag{2}$$

Consequently, the transfer function between current and acceleration is

$$G(s)=\frac{A(s)}{I(s)}=\frac{sV(s)}{I(s)}=sH(s)\tag{3}$$


Difference between diversity gain and coding gain

MIMO refers to any system with multiple antennas, where signals are transmitted over the same bandwidth and at the same time.

Space-time coding (block or trellis, or a different one) is a specific way to transmit the data symbols. In space-time block coding, symbols are repeated over time and space in a specific pattern, which creates (potentially) both transmit and space diversity. They don't offer an intrinsic coding gain because the information is not encoded.

Trellis codes, on the other hand, allow for a coding gain because not all symbol sequences are permissible. This is similar to non-MIMO trellis modulation, where previous transmitted symbols define a sub-set of the constellation that can be used for the following symbol. This is a form of encoding, with the main benefit of neither increasing the data rate nor the required bandwidth.


Choosing clock frequency and number of bits of an OFDM transceiver

The clock frequency $f_\mathrm{s}$ of the digital-to-analog converter (DAC) ultimately depends on the desired bitrate $R$ and some other parameters like guard interval length and number of bits per subcarrier. You can find the detailed formulas in a related question.

The speed of the IFFT block depends on $f_\mathrm{s}$ and the number of subcarriers $N_\mathrm{SC}$. It has to transfrom $f_\mathrm{S}/N_\mathrm{SC}$ sets of values each set consisting of $N_\mathrm{SC}$ complex numbers. What clock frequencies are involved in such an IFFT block depends on implementation details (mainly degree of parallelization) and cannot be answered in general.


Determining statistically the bandwidth of a signal from a spectrogram

A spectrogram is a tool to analyze non-stationnary signals, and hopefully identify or separate it onto simpler sources. Hence, parts of the spectrograms ought to be "different". Hence, from a statistical point of view, the 2D data is likely to be multimodal, and global statistical measures are likely to be useless.
If you now focus on specific "modes", like identifiable whistles in boxes, then you can retrieve 2D statistical spectral moments from those patches. These could be a 2D center of mass (in time and frequency) and the corresponding 2D spread (equivalent to a standard deviation) that would form an ellipsis around the center of mass. The vertical axis could provide an idea of the frequency span of this mode.
From my experience, this often requires to look at other statistical moments for more robustness.


Detecting a single frequency in noisy data

If your goal is to find only the frequency then its an easy task .
Magnitude = sqrt(real*real + img*img)
This number will always be positive. Now , find the maximum value of magnitude in the array.

also,
                   frequency = max magnitude*(Fs/N)

Fs = sampling Frequency(44100 in your case) , N = size of FFT
Now you have the frequency. But the numerical value depends upon the fft size that you are using i.e 512 size will result in accuracy of 86Hz (+- 86) , 1024 => 43Hz and so on. So , the frequency will be 189__ or 190__ something :)


How to plot BER Vs SNR of a received signal

From the plot it seems the transmitted QAM points are $(0.5, 0.5),(-0.5,0.5),(-0.5,-0.5),(0.5,-0.5)$, and if you can model your comm system as just AWGN channel, then your signal amplitude will just be $\sqrt{0.25+0.25} = \sqrt{0.5}$ and hence signal power will be $0.5$. In doing so, I have assumed all 4 QAM symbols equiprobable.

You can get your noise power by averaging all squared Euclidean distances of received points from their true QAM points.

From this you can get received SNR.


I want to actuate a motor based on a wav file, but it has too many samples per second

Well, if you are okay with a signal that follows the intensity instead of the envelope, which will be roughly the same, there is a simple trick to it.

Suppose you want one value every tenth of a second.  For each value, find the sample location at that time.  Multiply your signal by a suitable bell shaped curve centered at that point.  A VonHann window would work nicely.  Make sure the window interval is longer than several cycles of your sound signal.  Now, add the squares of the multiplied values over your windowed span.  Divide this sum by the number of samples in your window and take the square root.  This is the RMS (Root Mean Square, literally the square root of the mean of your squares) of your windowed signal.  Use this, or rescale it some, for your driving signal.

You will always get non-negative values.  Zero when the signal is zero, proportional to the amplitude for steady tones independent of frequency.

Here are the formulas:

$$ w[h] = \cos^2 \left( \frac{h}{2H} \pi \right) $$

$$ y[n] = \sqrt{ \frac{ \sum_{h=-H}^{H} ( w[h] x[n+h] )^2 }{ 2H+1 } } $$

The value of $H$ will determine how locally sensitive your measure is.  For H = 0, you get the absolute value of the signal.  Which suggests an alternative.  You could also simply use a weighted average of absolute values.  It would be comparable in results and more efficient to calculate.

The important thing is that $H$ be large enough to span several cycles of the underlying fundamental tone.

To make it even more efficient, you can skip a sample or two or three, as if your signal had a lower sampling rate.  But not too much, or results will suffer.


Complex Least Squares FIR filter design using Levinson' s algorithm

From this term in D1:

exp(-j*om1(231:430) *2)

it looks like you're asking it for a filter with a group delay of 2 samples even though the number of taps is 61.  Try setting this to:

exp(-j*om1(231:430) * 61/2)

i.e. set a desired group delay of about half the filter length.


Capacity and BER as performance of transmission system

Channel Capacity and Bit/Symbol Error rate are two different key performance indicators (KPI) of a system.
BER/SER, as it's obvious from its name, is  the number of Bit/Symbol errors divided by the total number of transferred Bits/Symbols. and you can use channel coding to improve it so it's not a fixed KPI.
while Channel Capacity is maximum rate that the system can achieve and it's a fixed KPI.
so it make sense to study the channel capacity of your system. by studying channel capacity you find maximum rate that can be achieved and by improving BER/SER you can improve rate of system to get close to channel capacity.


How do I identify a notch filter given the transfer function?

Remember that a notch filter must have transfer function zeros on the unit circle $z=e^{j\omega}$, i.e., there are certain frequencies that are completely attenuated. But note that zeros on the unit circle are not sufficient to characterize a filter as a notch filter (even though any frequency with perfect attenuation is called a notch frequency). E.g., most frequency-selective filters will have most, if not all, of their zeros on the unit circle (to get a good stop band attenuation). However, they are no notch filters. So for IIR filters you also have to consider the pole locations. And for notch filters, the pole angles are identical to or at least very close to the angles of the zeros. So for each zero on the unit circle you get a pole at the same or a similar angle, just with a smaller radius, because for causal and stable filters the poles need to be inside the unit circle.

That's why the filter in your example qualifies as a notch filter. All its zeros are on the unit circle and its poles are at the same angle as the zeros. The notch frequencies are simply the angles at which the zeros occur:

$$z_{01}=e^{j\omega_{01}}=e^{j\pi /4},\qquad
z_{02}=e^{j\omega_{02}}=-e^{j\pi/4}=e^{-j3\pi /4}$$


Hint:

Since $K(s) = \frac{1}{s-2+G(s)}$, then the roots of $s-2+G(s)$ must lie on the left side of the complex plane. Let's suppose that $G(s)$ is of the form $as + b$. Can you solve it then?


What is obtained from the cross correlation plot?

I may answer things that you already know,but i will do my best anyway (i had similar problems once..).
When you do the Correlation of the two sequences you simply shift the one relatively to the other (as you may already know) so each time you shift it you get lets say one "dot" for your plot over there,one single point.This points in a row construct the plot you are getting there.Now,in order to get each one of them you perform one operation,lets say better one calculation which is:the number of equal/same bits minus the number of non equal/same bits,after you perform that you shift again and calculate again....Now when it come to the point that two identical sequences are perfectly aligned,you can understand that the number of "non equal bits" will be zero (cause apparently will be matching perfectly and the result will be a massive number (compared to the rest) that will be YOUR CORRELATION PEAK,meaning that you have reached correlation.As i see there that is something similar to BOC for satellite coms that i have been using in the past,is that case that time   you need to shift the waveform in order to reach to the peak is the time that you need to uses to calculate the distance in your GPS system!
I hope that helps...


Your understanding is mistaken. LTI filters are applied via convolution in time domain and multiplication in frequency domain.

The confusion possibly happens about designing the filter vs applying it.

The windowed sinc filter impulse response is given by

$$ h[n] = w[n] h_i[n] ,$$

where $h_i[n]$ is the ideal (infinite length) sinc filter which corresponds to a brickwall lowpass type in your case. To approximate and taper that filter, you apply a window $w[n]$ of length $L$ to it by multiplication in time domain and its effect in the frequency domain is a convolution by the well known Fourier modulation theorem:

$$H(e^{j\omega}) =  W(e^{j\omega}) \star H_i(e^{j\omega}) $$

However, when the filter is applied on an input, the output equation in the time domain is

$$y[n] = h[n] \star x[n] ,$$

and the corresponding frequency domain equation is

$$Y(e^{j\omega}) = H(e^{j\omega}) X(e^{j\omega}).$$

So, rewrite your code base don this to get an expected result.


Finding the auto-correlation sequence $r_{xx}[k]$ for an AR(2) process

We have
\begin{equation}
x[n] = v[n] + 0.75x[n-1]-0.25x[n-2]
\end{equation}
You started correct, you get
\begin{equation}
 \begin{split}
  r_{xx}[k] = E\{(v[n] + 0.75x[n-1]-0.25x[n-2])(v[n-k] + 0.75x[n-k-1]-0.25x[n-k-2])\} \\
 \end{split}
\end{equation}
You get 9 terms
\begin{equation}
 \begin{split}
  r_{xx}[k] &= Ev[n]v[n-k] +0.75E x[n-k-1]v[n] - 0.25E x[n-k-2]v[n]\\
   &+0.75 Ex[n-1]v[n-k] + 0.75^2Ex[n-1]x[n-k-1]+0.75(0.25)Ex[n-1]x[n-k-2]\\
   &-0.25Ex[n-2]v[n-k] -0.75(0.25)Ex[n-2]x[n-k-1] + 0.25^2Ex[n-2]x[n-k-2]
 \end{split}
 \tag{1}
\end{equation}
Since the noise is white then
\begin{equation}
 Ev[n]v[n-k]  = \delta(k)\sigma_v^2 = \delta(k)
\end{equation}
Since $k\geq 0$, write down $x[n-k-1]$ and $x[n-k-2]$ do not contain lags that align with $v[n]$, so
\begin{align}
 E x[n-k-1]v[n] &= 0\\
 E x[n-k-2]v[n] &= 0
\end{align}
On the other hand, $Ex[n-1]v[n-k]$ and $Ex[n-2]v[n-k]$, will agree for $k = 1$ and $k=2$ according to
\begin{align}
x[n-1]v[n-k] &= \overbrace{Ev[n-1]v[n-k]}^{\delta(k-1)} + 0.75\overbrace{ Ex[n-2]v[n-k] }^{\delta(k-2)}- 0.25\overbrace{Ex[n-3]v[n-k]}^{0} \\
x[n-2]v[n-k] &= \underbrace{Ev[n-2]v[n-k]}_{\delta(k-2)} + \underbrace{E(0.75x[n-3]-0.25x[n-4])v[n-k]}_{0}
\end{align}
So
\begin{align}
 x[n-1]v[n-k] &= \delta(k-1) + 0.75\delta(k-2) \\
 x[n-2]v[n-k] &= \delta(k-2)
\end{align}
The other four terms follow from the AR definition, i.e.
\begin{align}
 E x[n-1]x[n-k-1] &= r_{xx}[n-1-n+k+1]= r[k] \\
 E x[n-1]x[n-k-2] &= r_{xx}[n-1-n+k+2]= r[k+1] \\
 E x[n-2]x[n-k-1] &= r_{xx}[n-2-n+k+1]= r[k-1] \\
 E x[n-2]x[n-k-2] &= r_{xx}[n-2-n+k+2]= r[k]
\end{align}
Replacing all results in equation $(1)$, we get
\begin{equation}
\begin{split}
  r[k] &= \delta(k) +0.75\delta(k-1) + 0.75^2\delta(k-1) -0.25\delta(k-2)\\
   &+0.75^2 r[k] +0.75(0.25)r[k+1]-0.75(0.25)r[k-1] + 0.25^2 r[k]
\end{split}
\end{equation}
which suggests that
\begin{equation}
 0.375r[k] = \delta(k) + 1.3125\delta(k-1) - 0.25\delta(k-2) +0.1875(r[k+1] - r[k-1])
\end{equation}
For $k = 0$, we get
\begin{equation}
 0.375r[0] = \delta(0) + 1.3125\delta(-1) - 0.25\delta(-2) +0.1875(r[+1] - r[-1])
\end{equation}
we know that $\delta(0) = 1$ and $\delta(-1) = \delta(-2) = 0$. Also $r[1] = r[-1]$ due to the fact that $r[k]$ is even
\begin{equation}
 r[0] =\frac{1}{0.375} \delta(0) = \frac{8}{3}
\end{equation}
For $k = 1$, we get
\begin{equation}
 0.375r[1] = 1.3125 + 0.1875(r[2] - r[0])
\end{equation}
which is (again because $r[k]$ is an even function)
\begin{equation}
 0.375r[1] = 1.3125 + 0.1875(r[2] - \frac{8}{3}) \tag{*}
\end{equation}
and at $k = -1$
\begin{equation}
 0.375r[-1] = 0.1875(r[0] - r[-2])
\end{equation}
which is (again because $r[k]$ is an even function)
\begin{equation}
 0.375r[1] = 0.1875(\frac{8}{3} - r[2])\tag{**}
\end{equation}
Now solve the system of two equations in two unknowns $(*),(**)$ to get $r[1],r[2]$, i.e.
\begin{equation}
 \begin{bmatrix}
    0.3750 &  -0.1875\\
    0.3750  &  0.1875
 \end{bmatrix}
 \begin{bmatrix}
  r[1]\\
  r[2]
 \end{bmatrix}
 =
 \begin{bmatrix}
    0.8125\\
    0.5000
 \end{bmatrix}
\end{equation}


Enlarging a picture by adding pixels in matlab

Yes, you can use imread to get the image values, interp2 for the interpolation, and imwrite to make a new file.  Pay attention to data format and scaling.  Read the docs for these 3 functions and you should be ok


First, observe that it's sufficient to show the property for a first-order all-pass filter, because any higher order all-pass filter can be written as a cascade of first-order all-pass sections. And if

$$|G_i(z)|<1\quad\text{for}\quad |z|>1$$

then also

$$\left|\prod_iG_i(z)\right|<1\quad\text{for}\quad |z|>1$$

must be satisfied.

Now derive an expression for the squared magnitude of the first-order all-pass transfer function $G_1(z)$:

$$|G_1(z)|^2=\frac{1-2\text{Re}\{\alpha^*z\}+|\alpha|^2|z|^2}{|z|^2-2\text{Re}\{\alpha^*z\}+|\alpha|^2}\tag{1}$$

Using $(1)$, the inequality $|G_1(z)|^2<1$ can be written as

$$\begin{align}1-2\text{Re}\{\alpha^*z\}+|\alpha|^2|z|^2&<|z|^2-2\text{Re}\{\alpha^*z\}+|\alpha|^2\\
1+|\alpha|^2|z|^2&<|z|^2+|\alpha|^2\\0&<(|z|^2-1)(1-|\alpha|^2)\tag{2}\end{align}$$

Note that if the all-pass filter is stable, we must have $|\alpha|<1$, because $\alpha$ is the pole location, and a causal and stable filter must have all its poles inside the unit circle of the complex plane. Consequently, the inequality $(2)$ is satisfied for $|z|>1$ because both factors on the right-hand side are positive.


Which is correct representation for function of variable?

Based on your new explanations, I think the answer is just
$$
E_{u_i}
$$
given that what you're really trying to find a shorthand for is the expression inside the summation sign in (5) of the screenshot. That is:

$$
\min_{U} \left\{ \lambda \sum_{i=1}^N \int_\Omega r_i(\chi, p_i(\chi)) u_i(\chi) d\chi \right\} = \min_{U} \left\{ \lambda \sum_{i=1}^N E_{u_i} \right\}
$$

so that
$$
E_{u_i} = \int_\Omega r_i(\chi, p_i(\chi)) u_i(\chi) d\chi
$$

Original answer below for reference

As per @Fat32's comment, any of:

$E(\sigma, z, p_i)$
$E_{\sigma,p_i}(z)$
$E(z; \sigma,p_i)$

is acceptable.

But you need to tighten up your notation. What is $x$ vs $z$? How does $\Omega$ relate to them? What is $y$?


Suppressing the lower sideband decreases the bandwidth by a factor of $2$ compared to BPSK. Since we're not talking about single-sideband (SSB) modulation, but vestigial sideband (VSB) modulation, we don't gain a factor of $2$, but a factor that is slightly smaller than two, because a part of one of the sidebands is retained.

Note that both BPSK and 2-VSB carry the information only in the (bi-polar) amplitude of the signal, i.e., only in the in-phase component, not in the quadrature component. The complex baseband signal for 2-VSB only occurs because the larger part of one of the sidebands is filtered out. This, however, is not necessarily done by complex baseband processing (using the Hilbert transform), but it can be done by applying a band pass filter after modulating the baseband signal. The latter option only uses real-valued processing.


How to convert magnitude to dB in matlab


Goertzel algorithm produces incorrect phase

The short answer is that when Goertzel algorithm is used for omega, which doesn't correspond to of the finite set of frequencies computed by DFT (which correspond to integer k), then an additional correction factor should be applied. I've got an answer from Peter K after describing my question more clearly here: What's wrong with my Goertzel algorithm implementation?


Pink ($1/f$) pseudo-random noise generation

Linear Filtering

The first approach in Peter's answer (i.e. filtering white noise) is a very straightforward approach. In Spectral Audio Signal Processing, JOS gives a low-order filter that can be used to produce a decent approximation, along with an analysis of how well the resulting power spectral density matches the ideal.  Linear filtering will always yield an approximation, but that may not matter in practice.  To paraphrase JOS:

  There is no exact (rational, finite-order) filter which can produce
  pink noise from white noise. This is because the ideal amplitude
  response of the filter must be proportional to the irrational function
  $1/\sqrt{f}$, where $f$ denotes frequency in Hz. However, it is easy
  enough to generate pink noise to any desired degree of approximation,
  including perceptually exact.

The coefficients of the filter he gives are as follows:

B = [0.049922035, -0.095993537, 0.050612699, -0.004408786];
A = [1, -2.494956002, 2.017265875, -0.522189400];

They're formatted as parameters to the MATLAB filter function, so for the sake of clarity, they correspond to the following transfer function:

$$
H(z) = { .041 - .096z^{-1} + .051z^{-2} - .004z^{-3} \over{}  1 -2.495z^{-1} + 2.017z^{-2} - .522z^{-3}}
$$

Obviously, it's better to use the full precision of the coefficients in practice. Here's a link to what the pink noise generated using that filter sounds like:

For fixed-point implementation, since it's usually more convenient to work with coefficients in the range [-1,1), some reworking of the transfer function will be in order.  Generally, the recommendation is to break things up into second-order sections, but part of the reason for that (as opposed to using first-order sections) is for the convenience of working with real coefficients when the roots are complex.  For this particular filter, all the roots are real, and combining then into second-order sections would probably still yield some denominator coefficients >1, so three first-order sections is a reasonable choice, as follows:

$$
H(z) = { 1 - b_1z^{-1}\over{1 - a_1z^{-1}} } \space { 1-b_2z^{-1} \over{1-a_2z^{-1} } } \space
{1 - b_3z^{-1} \over{ 1 - a_3z^{-1} } }
$$

where

$$
b_1 = 0.98223157, \space b_2 = 0.83265661, \space b_3 = 0.10798089
$$
$$
a_1 = 0.99516897, \space a_2 = 0.94384177, \space a_3 = 0.55594526
$$

Some judicious choice of sequencing for those sections, combined with some choice of gain factors for each section will be required to prevent overflow.  I haven't tried any of the other filters given in the link in Peter's answer, but similar considerations would probably apply.

White Noise

Obviously, the filtering approach requires a source of uniform random numbers in the first place.  If a library routine isn't available for a given platform, one of the simplest approaches is to use a linear congruential generator.  One example of an efficient fixed-point implementation is given by TI in Random Number Generation on a TMS320C5x (pdf).  A detailed theoretical discussion of various other methods can be found in Random Number Generation and Monte Carlo Methods by James Gentle.

Resources

Several sources based on following links in Peter's answer are worth highlighting.

The first filter-based chunk of code references Introduction to Signal Processing by Orfanidis.  The full text is available at that link, and [in Appendix B] it has coverage of both pink and white noise generation.  As the comment mentions, Orfanidis mostly covers the Voss algorithm.
The Spectrum Produced by the Voss-McCartney Pink Noise Generator.  Way down near the bottom of the page, after extensive discussion of variants of the Voss algorithm, this link is referenced in giant pink letters.  It's much easier reading than some of the preceding ASCII diagrams.
A Bibliograph on 1/f Noise by Wentian Li.  This is referenced both in Peter's source and by JOS.  It has a dizzying number of references on 1/f noise in general, dating all the way back to 1918.


Calculation of gain at certain frequency in case of nth order IIR filter

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.


HINT: The error energy can be written as
$$\begin{align}\sum_k\big|d[k]\big|^2&=\sum_k\left|x[k]-\frac{1}{2\pi}\int_{-W}^W\sum_nx[n]e^{-jn\omega}e^{jk\omega}d\omega\right|^2\\&=\sum_k\left|x[k]-\sum_nx[n]\frac{1}{2\pi}\int_{-W}^We^{j(k-n)\omega}d\omega\right|^2\tag{1}\end{align}$$
Now compute the integral
$$I(k-n,W)=\frac{1}{2\pi}\int_{-W}^We^{j(k-n)\omega}d\omega\tag{2}$$
and show that $\lim_{W\to\pi}I(k-n,W)=\delta[k-n]$. This reduces the second sum in $(1)$ to a single element and the result follows.


We can compute the marginal PDF of $X$ as follows:
\begin{align}
 f_X(x)=\int_\delta f_{X|D}(x|\delta)f_D(\delta)d\delta = f_0(x) \int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta
\end{align}

We can express the previous equation in terms of $\lambda(x)$ by dividing both sides by $f_0(x)$ and applying the natural logarithm.

\begin{align}
 \lambda(x) := \log \frac{f_X(x)}{f_0(x)} = \log \left[ \int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta\right]
\end{align}

Now, if we derive w.r.t the parameter $x$, we obtain the following equality:

\begin{align}
 \frac{\partial}{\partial x} \lambda (x) = \frac{\int_\delta \alpha \delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}{\int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}= \alpha \frac{\int_\delta \delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}{\int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}
\end{align}

With this equality, we are now able to express the MMSE estimator of $D$ in terms of $\lambda(x)$:

\begin{align}
 \mathbb{E}[D|X=x]&=\int_{\delta}\delta f_{D|X}(\delta|x)d\delta=\int_{\delta}\delta\frac{f_{X|D}(x|\delta)f_D (\delta)}{f_X(x)}d\delta=\frac{\int_\delta \delta f_0(x)e^{\alpha \delta x - \psi(\delta)}f_D(\delta) d\delta}{f_X(x)} = \\
 &= \frac{\int_\delta \delta f_0(x)e^{\alpha \delta x - \psi(\delta)}f_D(\delta) d\delta}{\int_\delta f_0(x) e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}=\frac{1}{\alpha}\frac{\partial}{\partial x} \lambda (x)
\end{align}


Short Answer
In order to do a complete Zero-IF transmit and receive without inversion, the operations must be done as follows (or otherwise have the signs equally flipped on Q for both transmit and receive):
Zero-IF (Direct Conversion) Transmit:
$$I_{RF}(t) = I_s(t)\cos(\omega_{LO}t) - Q_s(t)\sin(\omega_{LO}t)$$
Zero-IF Receive:
$$I(t) = \text{LPF}\{I_{RF}(t)\cos(\omega_{LO}t)\}$$
$$Q(t) = \text{LPF}\{-I_{RF}(t)\sin(\omega_{LO}t)\}$$
Yes spectral inversion can easily occur by having a sign change on Q on one side and not the other, showing that the way to correct it is to simply change the sign of the Q output.
Where:
$I_{RF}(t)$: Transmitter RF signal
$I_s(t)$: Baseband waveform I input at transmitter
$Q_s(t)$: Baseband waveform Q input at transmitter
$\cos(\omega_{LO}t)$: I LO
$\sin(\omega_{LO}t)$: Q LO
$\text{LPF}\{\}$: represents a low pass filtering operation.
Detailed Explanations
Consider the plots below that demonstrate exactly how this occurs:
For context first consider the typical Zero-IF Receiver architecture where we have a real signal that gets multiplied with a complex Local Oscilator, as implemented with two multipliers, one for the real or "in-phase" component as "I" and the other for the imaginary or "quadrature" component as "Q".  In the diagram the LO is shown to be fed into a quadature splitter (indicated by 90°) which produces two sinusoidal outputs in quadrature. The two outputs form a single complex source as real and imaginary components. One is designated to be "cosine" while the other is "sine" with "cosine" typically representing the real component. This is consistent with the Euler relationships listed further below.

Below shows the spectrums during this downconversion process. Note the Euler identities in the creation of the complex LO:
$$e^{+j\omega t} = \cos(\omega t) + j\sin(\omega t)$$
$$e^{-j\omega t} = \cos(\omega t) - j\sin(\omega t)$$
Multiplication in the time domain is convolution in the frequency domain, so the result of multiplying the real input with the complex LO (which is a single sided impulse in frequency) is to just shift the spectrum accordingly. So when we multiply by $e^{-j\omega t}$ the spectrum slides to the left as a result (and then low pass filter to get the resulting baseband signal), and when we multiply by $e^{+j\omega t}$, the spectrum slides to the right. Notice how the first case results in an inverted spectrum.

To reverse this inversion, we could simply change the sign of the sine (Q) output of the LO resulting in a negative frequency LO instead of a positive frequency. Alternatively and equivalently change the sign of the Q output of the receiver. To see this simply, consider a single baseband positive tone which in the time domain can be represented as a phasor rotating counter-clockwise on a complex IQ plane. If you negate Q (the imaginary component of the phasor), it would then rotate clock-wise instead, representing a single baseband negative tone of the same frequency.

The OP had additional questions specific to the transmitter side (direct-conversion transmitters for which I add the details below:
To explain the graphic we see how the transmitter structure can be viewed as the real component of a full complex multiplication; starting with the baseband Tx waveform on the top spectrum, the complex LO as $e^{+j\omega t}$ specifically to shift the spectrum to the right (for no inversion), resulting in what would be a complex RF signal on the third spectrum (which would require a real and imaginary datapath). By taking the real part of this, we get the conjugate symmetric real spectrum as on the bottom spectrum that can be fed to a single antenna.

To realize this, consider the complex baseband as $I_s+jQ_s$ multiplied with the complex LO as $e^{+j\omega_{LO}t}= \cos(\omega_{LO}t)+j\sin(\omega_{LO}t)$. The result of a full complex multiplication (requiring 4 real multipliers) would be:
$$(I_s(t)+jQ_s(t))\bigl(\cos(\omega_{LO}t)+j\sin(\omega_{LO}t)\bigr) $$
$$ = \bigl(I_s(t)\cos(\omega_{LO}t) - Q_s(t)\sin(\omega_{LO}t)\bigr)+j\bigl(I_s(t)\sin(\omega_{LO}t)+Q_s(t)\cos(\omega_{LO}t)\bigr)$$
With the real part representing the transmitted RF signal specifically as:
$$I_{RF}(t) = I_s(t)\cos(\omega_{LO}t) - Q_s(t)\sin(\omega_{LO}t)$$
Notice the minus sign on the Q LO input as required to not have an inverted spectrum. If one were to make the mistake in using $+\sin(\omega_{LO}t)$ for Q then a spectral inversion would result! This could very well be a source of confusion since $e^{+j\omega t} = \cos(\omega t)+j\sin(\omega t)$, but one needs to review carefully the operation above when taking the real of the full complex multiplication to see why the negative sign is required to actually realize the positive complex LO!
Similarly with the Rx, in order to not introduce and inversion, we must multiply the received real RF signal by $e^{-j\omega_{LO}t}$ as shown earlier in the first spectrum plot on the left hand side. But in this case the RF signal is completely real, so we can show as in the product done out for the transmitter but in this case use $I_{RF}$ as the signal multiplied with the complex LO as $e^{-j\omega_{LO}t}= \cos(\omega_{LO}t)-j\sin(\omega_{LO}t)$ to get the complex baseband output (prior to low pass filtering). This results in:
$$(I_{RF}(t))\bigl(\cos(\omega_{LO}t)-j\sin(\omega_{LO}t)\bigr) $$
$$ = I_{RF}(t)\cos(\omega_{LO}t) -jI_{RF}(t)\sin(\omega_{LO}t)$$
With the final zero-IF receiver output as:
$$I(t) = \text{LPF}\{I_{RF}(t)\cos(\omega_{LO}t)\}$$
$$Q(t) = \text{LPF}\{-I_{RF}(t)\sin(\omega_{LO}t)\}$$


Forward error correction with limited possible errors

Reed-Solomon codes are perfectly matched to a symmetric channel for which the probability that a given symbol $X$ is changed to a symbol $Y \neq X$ is the same for all choices of $Y \neq X$. That is, no particular symbol error is more likely than any other symbol error.

Thus, knowing that the channel is not symmetric and that certain symbol errors are impossible doesn't help in the decoding process at all. This knowledge might be usable in the post-decoder processing. If the decoder has corrected $Y$ into $X$ and $Y \notin S(X)$, then that is a dead giveaway that something is awry, right? So, the post-decoding processor can reject such purported "correct decoding" and ask for a retransmission or blank out that codeword to produce a blank area on the screen instead of gobbledygook.

Unfortunately, the most likely problem with Reed-Solomon decoders is not that they find an incorrect codeword, but rather that they fail to decode because the received word is not at Hamming distance $t$ or less from any codeword.

I don't think you are going to find any way of increasing the rate of the code based on your knowledge of the channel's asymmetry.


Multi-Stage digital FIR filter vs Single stage FIR filter. Which is better?

First, two cascaded 24th-order (25-tap) filter will yield a 48th-order filter which is not identical to a 49th-order (50 tap). However, as I assume that this isn't really the answer you are looking for, but rather the difference between multi-stage and single stage.

Assuming infinite precision arithmetic they are the same. In practice, multi-stage filters have lower coefficient sensitivity, so using two (or more) filters would (statistically) allow shorter coefficients word lengths. However, multi-stage filters also require intermediate quantization, so for a complete comparison one should also consider that. And, as obvious from the first paragraph, multi-stage require more multiplications for a given filter order. I would say that unless you run into problems with coefficient quantization you should go with a single stage.

Another use case where multi-stage is relevant is when implementing interpolation and decimation. By using more than one stage some of the filtering can be performed at a lower sample rate and the total arithmetic complexity is reduced.


What is different between gradient and imgradient functions in MATLAB

Question 1:
gradient returns both estimated horizontal $\left(\frac{\partial I}{\partial x}\right)$ and vertical components $\left(\frac{\partial I}{\partial y}\right)$ of the gradient vector, which in your example code are Gx and Gy. To do so, gradient computes for each point $X(i)$ the half difference of $X(i-1)$ and $X(i+1)$, which is the most basic gradient estimation method (it's called the central difference method).
Then, computing

g = sqrt(Gx.^2+Gy.^2);

gives you the gradient magnitude at each point.

imgradient differs since it directly outputs the gradient's magnitude (g) and direction (which you suppressed using ~) at each point. In addition, it computes the gradient using a two-dimensional method applying a sobel kernel to each point and its neighbors. Finally, this function allows you to specify which kernel to use to compute the numerical gradient if you don't want to use the Sobel estimator. You can for example choose to estimate the gradient with the central difference method as in gradient.

Question 2 :
Sobel method is just the simplest contour detection method providing good results for a low computation cost. But there are plenty of others, that are used depending on the application I guess.
Compared to the central difference method, it achieves better noise rejection using a different kernel with coefficients designed to do so. See as an example the kernel used for horizontal gradient estimation :
$$
\begin{matrix}
-1 & 0 & 1\\
-2 & 0 & 2\\
-1 & 0 & 1
\end{matrix}
$$
The considered point's row is emphasized (factor 2) compared to it's neighbors (factor 1). While the central difference method does not make that difference and weights all points with the same coefficient :
$$
\begin{matrix}
-1 & 0 & 1\\
-1 & 0 & 1\\
-1 & 0 & 1
\end{matrix}
$$


Creating odd harmonics without the fundamental

Im looking for a mathematical function which can recreate any fundamentals odd harmonics, in a multitone input...

(emphasis by me)

What you ask for is not possible. If you demand multiple tones and complex tones as they appear in music, then your problem leads to a conflict. A complex tone comes with overtones. Even if you have a process that adds odd harmonics to each sinusoidal component, then the even overtones of the fundamental frequency will contribute even harmonics. That is because an odd number multiplied with an even number is again even.

So circumvent that, you would have to identify the fundamental frequency of complex tones first, leading you directly to the problem of polyphonic pitch detection which remains mostly unsolved.

Maybe you should explain what you're really trying to accomplish, as in what the application is that you have in mind. With that information it's quite likely that we can come up with a relaxed set of requirements that can actually be implemented.


Estimate mean amplitude in frequency range using fft

The Discrete Fourier Transform is a summation :

$$ X_k =\sum_{n=0}^{N-1}x_ne^{-j2\pi kn/N} $$

Which means, the more data you have, the bigger the output mangitude. What you can do is normalizing your output by dividing the magnitude of the result by the number of samples in your signal.

I would invite you to check this question for more details.

On a side note, you may want to consider zero padding for better peak detection accuracy. Also, if you intend to filter out the noise by suppressing the peaks that you consider noise and then do an inverse FFT, you will gain to read this other question

We could most likely help you more if you gave more details about your problem as I suggested in my comment.

Hope that helps.


Change the melody of human speech using FFT and polynomial interpolation

The low quefrency MFCC coefficients might be a place to start.
I read a few articles by Stylianou from a few years back where he did funky stuff to peoples voice (gender and age modification).


Distorted output of delay VST effect in C

I am not familiar with signal processing for audio, but I would like give answer which contain common programming suggestions. Sorry, I really do not know your programming level and may be this suggestions will be look obvious and ridiculous, but I hope it will be useful.

  I'm getting some odd results with some code I'm working on.

When you get odd results usual this caused by one of three reasons:

incorrect implementation of correct logic;
correct implementation of incorrect logic;
incorrect implementation of incorrect logic.

Some implementation errors such as: using of uninitialized variables, using of wrong pointers, accessing to memory out of bounds can be easy found when using some special tools. One kind of tools is static code analyser, there are a large number of varied of static code analysers, but I would suggest you use cppcheck. Other kind of tools is (memory) profilers. I would suggest you use memcheck tool from valgrind suite. Very often this tools helps find very tricky errors.

To find logical error you need test code with different input data and conditionals. It is possible to do it manually time to time, but more easy to do it automatically by calling small code snippets. This code snippets can be designed as simple console application or as test case using some testing framework. Testing framework is not necessary for writing test, but it need for manage tests. I suggest to use the C++ Google testing framework.

Now I will try to implement my own suggestions.

First of all to allow a piece of existed code will be compiled without errors I added header file and put into it all required types and variables:

#include <cstdint>

using VstInt32 = int32_t;
constexpr int delay_samples             {1};
constexpr int CIRCULAR_BUFFER_LENGTH    {4};
constexpr float feedback             {0.0f};
constexpr float delay_level          {0.0f};

struct CircularBuffer
{
    CircularBuffer ()
    {
        for (int i = 0; i < CIRCULAR_BUFFER_LENGTH; ++i) {
            buffer [i] = 0.0f;
        }
        index = 0;
    };

    int index;
    float buffer [CIRCULAR_BUFFER_LENGTH];
};

class delay
{
    public:
        void processReplacing (float **, float **, VstInt32);

    private:
        CircularBuffer cBuf;
};

Shared code I placed into source (*.cpp) file:

#include "delay.h"

void delay::processReplacing (float ** inputs, float ** outputs, VstInt32 sampleFrames)
{
    float * in  = inputs  [0];
    float * out = outputs [0];
    float delayed_sample = 0.0f;
    int delayed_sample_index;
    for (int frame_idx = 0; frame_idx < sampleFrames; frame_idx++) {
        delayed_sample_index = (((int)(cBuf.index - delay_samples) + CIRCULAR_BUFFER_LENGTH) % (int)CIRCULAR_BUFFER_LENGTH);
        delayed_sample = cBuf.buffer [delayed_sample_index];
        // Write current sample to buffer
        cBuf.buffer [cBuf.index] = in [frame_idx] + (delayed_sample * feedback);
        // Calculate + write output sample
        out [frame_idx]          = in [frame_idx] + (delayed_sample * delay_level);
        // Increment buffer index
        cBuf.index = (cBuf.index + 1) % CIRCULAR_BUFFER_LENGTH;
    }
}

Now time to write tests. I am not familiar with logic of function which you try to implement, so I write one simple test just for example:

#include <gtest/gtest.h>
#include <iostream>
#include "delay.h"

class delay_test : public ::testing::Test
{
    public:
        void SetUp    () override;
        void TearDown () override;

        static constexpr int sampleFrames {8};
        float ** inputs;
        float ** outputs;
};

// This code will be automatically called before every test.
// This is feature of Google Tests Framework.
void delay_test::SetUp ()
{
    inputs      = new float * [1];
    outputs     = new float * [1];
    inputs  [0] = new float   [sampleFrames];
    outputs [0] = new float   [sampleFrames];
    for (int i = 0; i < sampleFrames; ++i) {
        inputs  [0][i] = 0.0f;
        outputs [0][i] = 0.0f;
    }
}

// This code will be automatically called before every test.
// This is feature of Google Tests Framework.
void delay_test::TearDown ()
{
    delete [] inputs  [0];
    delete [] outputs [0];
    delete [] inputs;
    delete [] outputs;
}

TEST_F (delay_test, response_on_1_Success)
{
    // Arrange.
    delay d;
    inputs [0][0] = 1.0f;

    // Act.
    d.processReplacing (inputs, outputs, sampleFrames);

    // Assert.
    for (int i = 0; i < sampleFrames; ++i) {
        if (1 == i) {
            EXPECT_FLOAT_EQ (1.0f, outputs [0][i]);
        }
        else {
            EXPECT_FLOAT_EQ (0.0f, outputs [0][i]);
        }
    }
}

Last step is writing build script. I use cmake, so simple build script has the following content:

cmake_minimum_required (VERSION 2.8)

project (delay)

set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -mtune=native -std=c++11 -Wall -Wextra")

set (
    delay_test_SOURCES
    ${PROJECT_SOURCE_DIR}/delay.cpp
        ${PROJECT_SOURCE_DIR}/delay_test.cpp
)

include_directories (
    ${PROJECT_SOURCE_DIR}
)

add_executable (
    delay_test
    ${delay_test_SOURCES}
)

target_link_libraries (
    delay_test
    gtest
    gtest_main
    pthread
)

Now it is possible to test code.

Let's check code with cppcheck:

$ cppcheck -I. --enable=all *.cpp
Checking delay.cpp...
[delay.cpp:7]: (style) The scope of the variable 'delayed_sample' can be     reduced.
[delay.cpp:8]: (style) The scope of the variable 'delayed_sample_index' can be reduced.
1/2 files checked 33% done
Checking delay_test.cpp...
2/2 files checked 100% done
Checking usage of global functions..
[delay_test.cpp:19]: (style) The function 'SetUp' is never used.
[delay_test.cpp:33]: (style) The function 'TearDown' is never used.
(information) Cppcheck cannot find all the include files (use --check-config for details)

The results are good, nothing suspicious is found.

Now let's launch tests and simultaneously memory profiler:

$ valgrind ./delay_test
==11407== Memcheck, a memory error detector
==11407== Copyright (C) 2002-2013, and GNU GPL'd, by Julian Seward et al.
==11407== Using Valgrind-3.10.0 and LibVEX; rerun with -h for copyright info
==11407== Command: ./delay_test
==11407==
Running main() from gtest_main.cc
[==========] Running 1 test from 1 test case.
[----------] Global test environment set-up.
[----------] 1 test from delay_test
[ RUN      ] delay_test.response_on_1_Success
/home/gluttton/circ_buf/delay_test.cpp:57: Failure
Value of: outputs [0][i]
  Actual: 1
Expected: 0.0f
Which is: 0
/home/gluttton/circ_buf/delay_test.cpp:54: Failure
Value of: outputs [0][i]
Actual: 0
Expected: 1.0f
Which is: 1
[  FAILED  ] delay_test.response_on_1_Success (81 ms)
[----------] 1 test from delay_test (95 ms total)

[----------] Global test environment tear-down
[==========] 1 test from 1 test case ran. (144 ms total)
[  PASSED  ] 0 tests.
[  FAILED  ] 1 test, listed below:
[  FAILED  ] delay_test.response_on_1_Success

 1 FAILED TEST
==11407==
==11407== HEAP SUMMARY:
==11407==     in use at exit: 0 bytes in 0 blocks
==11407==   total heap usage: 306 allocs, 306 frees, 61,494 bytes allocated
==11407==
==11407== All heap blocks were freed -- no leaks are possible
==11407==
==11407== For counts of detected and suppressed errors, rerun with: -v
==11407== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)

Good news is that no errors was found by valgrind (we can realize it by valgrind output messages). But bad news is that test is fell, it means that code works not as I expect.

Conclusions: in my opinion your code implemented without visible implementation errors, but you need check logic of your code with various input data.


How to choose order and cut-off frequency for low-pass Butterworth filter?

an $n$th-order Butterworth filter in the $s$-domain has a magnitude function that is two straight lines connected with a soft corner at $\omega_0$ and -3 dB.  the straight line at the left is a flat 0 dB but the line on the right (assuming you're looking at log frequency and dB) has a slope of $-20n$ dB/decade drop.  That is the same as $-6n$ dB/octave drop, if you like your log frequency in units of octaves.  So you have to figure out if the drop is steep enough too whack the frequency components you want whacked.

$$ \Big| H_n(j \omega) \Big|^2 = \frac{1}{1 + \left(\frac{\omega}{\omega_0}\right)^{2n}} $$

$$\begin{align}
 20 \log_{10} \big| H_n(j \omega)\big| &=-10\log_{10}\left(1 + \left(\tfrac{\omega}{\omega_0}\right)^{2n}\right) \\
\\
&\approx  \begin{cases} 0 \qquad & \omega \ll \omega_0 \\ -10\log_{10}\left( \left(\tfrac{\omega}{\omega_0}\right)^{2n}\right) & \omega_0 \ll \omega \end{cases} \\
\\
&=  \begin{cases} 0 \qquad & \omega \ll \omega_0 \\ 20 n \log_{10}( \omega_0) - 20 n \log_{10}( \omega) & \omega_0 \ll \omega \end{cases} \\
\\
\end{align}$$


But if we see this new vector we can see that here no translation of the signal is done.

You have modified the data in the frequency domain, by multiplying by exp(−j2πkD/N). The translation thus occurs in the time domain. It would be indeed weird if multiplyig a vector by exp(−j2πkD/N) would shift this same vector!

The inverse DFT of [11,3-2j,-1,3+2j] is: [4, 4, 1, 2].

The inverse DFT of [11,-3+2j,-1,-3-2j] is: [1, 2, 4, 4].

The shift property is observed - the signal has been shifted by 2 (modulo 4).


(FFT) Is it good practice to ignore each bin's phase when characterizing a signal

It depends,

sometimes phase matters and sometimes it doesn’t.

There is a paper by Oppenheim.

  Oppenheim, Alan V., and Jae S. Lim. "The importance of phase in
  signals." Proceedings of the IEEE 69.5 (1981): 529-541.

where an image is reconstructed with just phase and with just magnitude. The phase only image is nearly the same while amplitude alone is a mess.

For something like a stationary random noise, the signal model explicitly notes that the phase is random and is ignored.

There are probably a lot of cases where people ignore phase because it’s easier. Phase unwrapping can be difficult.  Interpretation can be difficult.

For cross correlation phase is very important.

It's never good practice to do things blindly.


How does spatially white noise differ from diffuse noise

In a diffuse noise field the noise energy flows in all directions with equal probability. This does not mean that the noise at one point is uncorrelated with the noise measured at another point. On the other hand, for spatially white noise, also called incoherent noise, the noise signals at two different points in space are always uncorrelated.

An example for a diffuse noise field is noise in a car or in an office. An example for spatially white noise would be amplifier noise in a microphone array, where the noise in each channel is independent of the noise in other channels.


Appropriate sampling frequency for bpsk modulated signal

The sampling frequency depends not only on the carrier frequency and the symbol (baud) rate, but also on the pulse shape. If you transmit with square pulses, you'll need a very large bandwidth (in theory, infinite). If you transmit with sinc pulses, then the base-band bandwidth is half the pulse rate. Other pulses will fall somewhere in between.

To expand on Peter K.'s comment: you should look into using the complex envelope for simulation. The advantage is that, when you use it, the sampling frequency depends only on the symbol rate and pulse shape, not on the carrier frequency.


Why isn't the Hartley Transform more widely used?

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

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

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

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


Hilbert transform of symbolic function in Matlab

The real part, the Hilbert transform and the derivation are all operators that commute.  So you can move the derivative just close to the function.


minimum sampling rate for very short duration signals

Any signal of finite extent (duration) in one domain (time) has infinite extent (bandwidth) in the other (frequency).  The bandwidth of a short pulse above some noise floor can be extremely high, due to the pulse window's transform convolved with the carrier frequency.

Sharp edges or cutoffs imply very high bandwidths.  Thus in practice, most designed pulses are shaped closer to Gaussian, with much wider tails than the main pulse lobe (in both domains).

Even so, if you know the exact frequency of an unmodulated fragment of a pure sine wave, then in zero noise, 2 or 3 non-aliased points are theoretically enough to solve the equations for the sinewave's  amplitude and phase.


Yes indeed modern signal processing uses matrix methods of linear algebra (or Linear System Theory more correctly), in addition to the classical calculus and harmonic analysis.

Linear algebra (matrix) based DSP considers (1D) signals as vectors and transforms as matrices that map inputs to outputs. As CedronDawg also mentioned, in such a representation the DFT (discrete Fourier transform) is merely a matrix mapping:

$$ Y = \Phi X $$ where $X = [x_1,x_2,...,x_N]^T$ is an $N \times 1$ input vector, $Y = [y_1,y_2,...,y_N]^T$ is the output vector,  and $\Phi$ is the $N \times N$ DFT transform matrix.

For a brief treatment of matrix methods in DSP, have a look at the books:

Multiresolution Signal Decomposition_AKANSU
Statistical Digital Signal Processing_HAYES


As you've mentioned, you can simply add the impulse responses, and, consequently, the frequency responses of LTI systems connected in parallel. In general, frequency responses are complex-valued, so you can't just add their magnitudes. However, if their phase responses are identical, you can add their magnitudes. A special case are two linear-phase FIR filters of the same length, which have the same (linear) phase response, so their magnitudes can be added.

If both filters approximate a magnitude of $1$ in their passbands, and a magnitude of zero in their stopbands, then the resulting filter will approximate a magnitude of $2$ where their passbands overlap, a magnitude of $1$ where the passband of one filter overlaps with the stopband of the other, and a magnitude of zero where both stopbands overlap. Note that FIR filters cannot have a constant magnitude (unless they have only $1$ filter tap), so all constant values are just approximated.

Now that you've figured out the resulting response when the outputs of the two filters are added, maybe you can come up with a creative way to connect the outputs of two linear phase FIR low pass filters (of the same length) with different cut-off frequencies to implement a bandpass filter.


$ {L}_{0} $ Pseudo Norm Minimization in Compressive Sensing

the solution for a sparse recovery problem is given by:
$$\text{min} ||x||_0$$
                    $$\text{s.t} \hspace{2mm}  y = Ax$$

The definition of $||x||_0$ is no. of non-zero entries in $x$. This is also called the sparsity of the vector.

i.e., we are asking for the sparsest solution $x$, that satisfies $y = Ax$.

Consider the simplest case where $x$ is 1-sparse but you don't know the location of that non-zero entry. In such a case, we have $\binom{N} {1}$ possibilities for a 1-sparse vector and to find the solution we have to examine the values of all the $\binom{N} {1}$ possibilities for finding the unique minimizer.
Similarly, if you are told that $x$ is $k$-sparse, you need to search $\binom{N} {k}$ possibilities.
i.e., the algorithm grows as $\binom{N}{k}$ with increase in $k$. Since $k$ is not known a priori, you have to check for all the $N$ possible values of $k$.

Hence, the complexity of the algorithm is $\displaystyle\sum_{i=1}^N \binom{N}{i}$. This is called as NP hard problem which is an acronym for Non-polynomial time complexity.


Digital RF signal vector resize algorithm

Probably not what you want, but using the convex hull of your points will preserve the minima and maxima.

R Code Below

#31780
T <- 1000
x <- rnorm(T,0,1)

ix <- seq(1,T)

ix_hpts <- chull(x = ix, y= x)

plot(ix_hpts, x[ix_hpts], type='l', col='red')
points(ix, x)


Well, you're not going to send and receive those symbols directly, you know?

What you'll need to convey the information over a channel is a waveform. Your receiver will then typically use the same waveform to correlate the signal. Say you are using root raised cosine pulses to transmit, then your receiver will employ a root raised cosine filter as a pulse matched filter to maximize the SNR.

There may be other things you can matched filter against. For instance, if you have an idea what your channel looks like, you might use a channel matched filter to maximize the SNR conveyed over your channel. Or maybe you are using spread spectrum techniques to transmit your information or things like code multiplex? You might want to use a code matched filter.

So, matched filtering can be used for a lot of things, the term alone doesn't tell it...


Evaluation of phase difference between two acquired signals (expected to be sine waves) with the same frequency

I really dislike the convention of notation in the question.  "$f$" is for frequency.  Functions having capital letters are for Fourier or Laplace transforms of the time-domain signal having the lower-case letter denoting it.

Let's say our reference signal is

$$ x[n] \triangleq A_x \cos(\omega n + \theta_x) $$

Somehow we can get the Hilbert transform of that signal

$$ \hat{x}[n] \triangleq \mathscr{H} \big \{ x[n] \big \} = A_x \sin(\omega n + \theta_x) $$

If you're generating the reference signal $x[n]$, you should also be able to synchronously generate $\hat{x}[n]$.   If you're measuring $x[n]$, you have to make a hilbert transformer to create $\hat{x}[n]$ (which will be a special FIR filter) and that will introduce a known delay.  You must delay all other signals by the same amount to get this relative phase measurement right.  The general relation of the discrete Hilbert transform is:

$$ \hat{x}[n] = \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m] $$

Of course, you will need to window that infinite sum to a finite sum and delay it to make it causal.  Note that all even $m$ terms in the above summation are zero.

$$\begin{align}
 \hat{x}[n] &= \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\
 &\approx \sum\limits_{m=-2M+1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\
 &= \sum\limits_{m=1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} (x[n-m]-x[n+m]) \\
 &= \sum\limits_{m=0}^{M-1} w[2m+1]\frac{2}{(2m+1) \pi} (x[n-2m-1]-x[n+2m+1]) \\
 &= \sum\limits_{m=0}^{M-1} h[m] \, (x[n-2m-1]-x[n+2m+1]) \\
\end{align}$$

where $w[n]=w[-n]$ is a decent window (like a Kaiser) of width $4M-1$ and the hilbert coefficients are

$$ h[n] \triangleq \frac{2\,w[2n+1]}{(2n+1) \pi} $$

Since this symmetrical, acausal FIR filter looks $2M-1$ samples into the "future", it must be delayed by $2M-1$ samples to be causal and realizable in a real-time context.  So you can not realize in real-time a true Hilbert transformer, but you can realize a very close approximation with the corresponding delay of $2M-1$ samples.  You must delay all other signals by the same amount to make the following calculations valid.  This means you will really be working with $\hat{x}[n-2M+1]$ and you have to explicitly delay $x[n]$ and $y[n]$ by the same amount.

So, while in real-time reality you will be working with $x[n-2M+1]$, $\hat{x}[n-2M+1]$, and $y[n-2M+1]$, for clarity, I am not showing the signals below as so delayed.

Now the measured signal that has the same frequency $\omega$ but different phase is:

$$ y[n] = A_y \cos(\omega n + \theta_y) $$

and the number we are looking for is $\theta_y - \theta_x$ .

Consider these two products:

$$\begin{align}
y[n] \cdot x[n] &= A_y A_x \cos(\omega n + \theta_y) \cos(\omega n + \theta_x) \\
   &= \frac{A_y A_x}{2} \big( \cos(2\omega n + \theta_y+\theta_x) + \cos(\theta_y - \theta_x)  \big) \\
\\
y[n] \cdot \hat{x}[n] &= A_y A_x \cos(\omega n + \theta_y) \sin(\omega n + \theta_x) \\
   &= \frac{A_y A_x}{2} \big( \sin(2\omega n + \theta_y+\theta_x) - \sin(\theta_y - \theta_x)  \big) \\
\end{align}$$

Like you have done before, you want to compute the mean of these two products:

$$\begin{align}
 \overline{y[n] \, x[n]} & = \frac{A_y A_x}{2} \cos(\theta_y - \theta_x) \\
\\
 \overline{y[n] \, \hat{x}[n]} & = -\frac{A_y A_x}{2} \sin(\theta_y - \theta_x) \\
\end{align}$$

We get the mean by use of any decent low-pass filtering.  We don't even need to worry that the DC gain = 1 as you normally do with the mean, since the scaling will come out in the wash.
Then, using complex expression:

$$\begin{align}
 e^{j(\theta_y - \theta_x)} &= \cos(\theta_y - \theta_x) + j \, \sin(\theta_y - \theta_x) \\
  &= \frac{2}{A_y A_x} \left( \overline{y[n] \, x[n]} - j \, \overline{y[n] \, \hat{x}[n]} \right) \\
\end{align}$$

and

$$\begin{align}
\theta_y - \theta_x &= \arg \bigg\{ \overline{y[n] \, x[n]} - j \,\overline{y[n] \, \hat{x}[n]} \bigg\} \\
 &= -\operatorname{atan2}\left(\overline{y[n] \, \hat{x}[n]}, \ \overline{y[n] \, x[n]}\right) \\
\end{align}$$

To do the four-quadrant arctangent:

$$\begin{align}
\operatorname{atan2}(y,x) &= \begin{cases}
\arctan(\frac y x) &\text{if } x > 0, \\
\arctan(\frac y x) + \pi &\text{if } x < 0 \text{ and } y \ge 0, \\
\arctan(\frac y x) - \pi &\text{if } x < 0 \text{ and } y < 0, \\
+\frac{\pi}{2} &\text{if } x = 0 \text{ and } y > 0, \\
-\frac{\pi}{2} &\text{if } x = 0 \text{ and } y < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\\
&= \begin{cases}
\arctan\left(\frac{y}{x}\right) &\text{if } x > 0, \\
\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y > 0, \\
-\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y < 0, \\
\arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\end{align}$$


Why does downsampling increase the phase angle of poles in autoregressive models?

Decimating a signal $x[n]$ by an integer factor $M$ results in the following $\mathcal{Z}$-transform:

$$Y(z)=\frac{1}{M}\sum_{k=0}^{M-1}X\left(e^{-j2\pi k/M}z^{1/M}\right)\tag{1}$$

If you assume that $x[n]$ was filtered by an anti-aliasing filter before decimation, the transfer function of the decimated filter can be approximated by the term for $k=0$ in (1):

$$Y(z)\approx\frac{1}{M}X(z^{1/M})\tag{2}$$

So if $z_{\infty}$ is a pole of $X(z)$ then $z_{\infty}^{M}$ is a pole of $Y(z)$. Writing the pole in polar coordinates $z_{\infty}=re^{j\phi}$ gives

$$z_{\infty}^M=r^Me^{jM\phi}$$

So the angle of the pole is increased by a factor $M$.


A snippet from the MATLAB documentation of the plot command:

  plot(Y) creates a 2-D line plot of the data in Y versus the index of each value...If Y is complex, then the plot function plots the imaginary part of Y versus the real part of Y, such that plot(Y) is equivalent to plot(real(Y), imag(Y)).

The circular plot you're seeing is corresponding to the real part on the x-axis and the imaginary part of the y-axis.

If you're doing plot(s), then MATLAB is plotting each column in s.


compute RMS from audio signal to get power
do AGC (automatic gain control)
perform "discrete differentiation" (the simplest is 1st order: $y[i] = x[i] - x[i-1]$)
if the value is greater than certain threshold, it means we have an onset. You have to determinate the threshold experimentally or use adaptive algorithm. Obviously you also need some kind of lowpass filter to filter the noise.

You may want to swap steps 2 and 3.

If you can use software libraries, check out aubio.

Possible duplicate of https://stackoverflow.com/questions/294468/note-onset-detection


Measuring walking speed with iPhone accelerometer

I don't have any idea about Accelerometer in iPhone, but converting Accelerometer data into distance is fairly difficult because the offset in double integration causes drift as already told. Check this interesting presentation sensor fusion for further details.
As far as I  know this can be avoided by using Digital Accelerometers( with high sampling rates), double integration has to be performed using electronic integrators with appropriate filters at every stage to remove the offset as well as the low frequency noise.
Also in that presentation he mentions the possibility of using more sophisticated algorithms like use of Kalman Filter for object
tracking or Pedometer Algorithm for pedometer design.

Hope this gives little more insight into your problem:)

PS: This does not accurately address your question still I thought I will share my knowledge, I want to comment but I don't have enough reputation !


The System Identification Toolbox contains facilities for

Linear system identification and determination of transfer function
Non-linear "black-box" models
ODE parameter estimation
Time series identification

One of the more promising substitutes for all of this functionality would be to use R.  Two packages that incorporate some of the time series functionality are the Time Series Utilities and Analysis package and the Flexible Modeling Environment for  Inverse Modelling, Sensitivity, Identifiability, Monte Carlo Analysis.  Also sharing these functions was the Time Series Analysis and Control package

An option for the transfer function calculation and manipulation (which attempts to be a substitute for the Matlab Control Systems toolbox, which overlaps with the System ID one significantly) is the Control Systems Library for Python, but the CSL doesn't have any functions for dealing with system identification, just the ability to manipulate transfer functions (thanks to JasonS for pointing this out).  The Ch Control System Toolkit allows one to model linear systems and determine the transfer function.

I'm sure there are other examples, but they may take some digging to uncover.


How to get general projection of rect function?

Since your expression
$$\begin{align*}P(\theta,t) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Pi(y)\Pi(x)\delta\left(x\cos\theta + y\sin\theta -t\right)dxdy\\
\end{align*}$$

uses a delta function with only 1 dimensional argument, in the 2 dimensional integration I assume you mean an infinite impulse line (an impulse "wall") and not a point impulse (an impulse "spike").

Let's just look at the general solution first, and then handle the degenrate case of $\theta = 0$.

I'll be using the following properties of the delta function:

$$\begin{align*}\delta(ax+b) &= \dfrac{1}{|a|}\delta\left(x + \dfrac{b}{a}\right)\\
\\
\int_{-\infty}^{\infty} f(x)\delta(x-a)dx &= f(a)\\
\end{align*}$$

So, I believe the solution is this:

$$\begin{align*}P(\theta,t) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Pi(y)\Pi(x)\delta\left(x\cos\theta + y\sin\theta -t\right)dxdy\\
\\
&= \int_{-\infty}^{\infty} \Pi(y) \left[\int_{-\infty}^{\infty} \Pi(x)\dfrac{1}{|\cos\theta|}\delta\left(x - \left[-\dfrac{\sin\theta}{\cos\theta}y +\dfrac{t}{\cos\theta}\right]\right)dx\right]dy\\
\\
&= \dfrac{1}{|\cos\theta|}\int_{-\infty}^{\infty} \Pi(y) \Pi\left(-\dfrac{\sin\theta}{\cos\theta}y +\dfrac{t}{\cos\theta}\right)dy\\
\\
&\mbox{and if the rectangles have any overlap for the given $t$ and $\theta$}\\
\\
&= \dfrac{1}{|\cos\theta|}\int_{\max\left(-\frac{1}{2},-\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}^{\min\left(\frac{1}{2},\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)} dy\\
\\
&= \dfrac{y}{|\cos\theta|}\biggr\rvert_{\max\left(-\frac{1}{2},-\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}^{\min\left(\frac{1}{2},\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}\\
\\
P(\theta,t) &= \dfrac{\min\left(\frac{1}{2},\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)-\max\left(-\frac{1}{2},-\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}{|\cos\theta|}
\end{align*}$$

So for the degenerate case of $\theta = 0$, the above expression runs into problems, but the integral for this special case is a bit easier:

$$\begin{align*}P(0,t) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Pi(y)\Pi(x)\delta\left(x\cos{0} + y\sin{0} -t\right)dxdy\\
\\
&= \int_{-\infty}^{\infty} \Pi(y) \left[\int_{-\infty}^{\infty} \Pi(x)\delta\left(x-t\right)dx \right]dy\\
\\
&= \Pi(t)\int_{-\infty}^{\infty} \Pi(y) dy\\
\\
&= \Pi(t)
\end{align*}$$

The one dimensional Fourier Transform of $P(0,t)$ is a simple table look-up to find that it is a $\mathrm{sinc}()$ function.

The one dimensional Fourier Transform of the more general $P(\theta, t)$ is obviously more complicated.  There are 4 cases there.


FIR filter design with matlab firls() function

The function firls() is meant to design filters with piecewise constant magnitude responses. So in practice you use only a few frequency points and the corresponding desired magnitude values, and the function computes a linear interpolation between the given frequency points. Of course, in theory your call to firls is correct, but I guess that the resulting system of linear equations becomes ill-conditioned.
I wrote a function lslevin.m, which can be used in the way you intended to use firls:

h = [0.8367 -1.537 1 -1.537 0.8367];
[H,w] = freqz(h,1,2048);    % that's more than enough frequency points
h2 = lslevin(5,w,H,ones(size(H)));
[H2,w] = freqz(h2,1,2048);
plot(w/pi,abs(H),w/pi,abs(H2),'r--')


Discrimiation between human speech and computer generated one

Many current computer speech synthesizers seem to produce less jitter or variability in inter glottal closure periodicity timing than most typical human speech.  (e.g. humans often produce subtle vibrato as well as tremolo, but synthesizers could change to do this as well.)


How to solve this Laplace integral for an averaging function in an iterative numerical (finite difference) model?

The variable of the weighted mean value zm would convey a more expressive power if denoted as zmean. The RHS integral in formula (14) of your reference is a "g(x)-weighted average of f(x) on the interval [0, t]", see http://www.sosmath.com/calculus/integ/integ04/integ04.html. In formula (14), f(s) is z2(s) and g(s) = sn. The integral in (14) is not a Laplace transform, but the extended g(x)-weighted mean value, and s is not a dual variable of the Laplace transform, but a dummy variable of integration.
A variable n in (14) is not an index but the exponent of the power function g(s) = sn, and t is also not indexed, but raised to a power n+1. From the fragment you presented in your question, and if they indeed use Newark's method, it does not follow that the authors seek an iterative solution for zm; rather, they use the weighted means as per formula (14) to find the solution for x, x_dotted.


Fourier Transform of a PSD and response of a PSD input

The Fourier Transform of stationary white noise exists: in fact "white" here refers to the characteristics of the Fourier Transform, in that the power spectral density will be constant over the entire spectrum, and also implies that each sample in the time domain is independent of all other samples.
The distribution of the magnitude of the samples in time is also usually defined, such as being "Gaussian" or "Uniform" depending on what distribution is actually used. For white noise samples the phase is also distributed uniformly  $\in [0, 2\pi)$. Due to the central limit theorem, the Fourier Transform of any distribution in time that is white will tend toward Gaussian in frequency. The resulting Fourier Transform will be Gaussian distributed in magnitude, uniformly distributed in phase, and each sample in the frequency domain is independent of all other samples.
The Power Spectral Density is computed from the Fourier Transform as the complex conjugate product, and not the other way around. The reason is once computed all phase information is lost so there is a many to one relationship in the computation.


Modulation Index of Conventional AM Signal

Modulation index for linear (continuous wave) amplitude modulation may be defined in different ways. But the one I found relevant to your definition is:

$$ \mu = \frac{ | \min \{ m(t) \} | } {A_c} $$

where $A_c$ is the carrier amplitude, $m(t)$ is the message signal and $\mu$ is the modulation index.

For simple demodulation of a conventional AM modulated signal $s(t)$ we look for $\mu < 1$ so that the envelope of $s(t)$ never becomes less than zero.

$$s(t) = [ A_c + m(t) ] \cos(2 \pi f_0 t) $$

Your signal is:
$$ s(t) = [ 100 + 100 \left( ~~\sin(2000 \pi t)+5\cos(4000 \pi t)~~ \right) ] \cos(2 \pi f_c t) $$

Your message signal is:
$$ m(t) = 100 \left( ~~\sin(2000 \pi t)+5\cos(4000 \pi t)~~ \right) $$

whose min value is $-600$ and from which you deduce the modulation index as:

$$ \mu = \frac{ | \min \{ m(t) \} | } {A_c} = \frac{ | -600 | } {100}  = 6 $$

I'm sorry to say but this is severely overmodulated...


Detecting Zero Crossing in Image Filtered by Laplacian Operator


How to distinguish between received quadratures

You would need to recover some timing information from the transmitted signal.

A common way is to transmit some sort of preamble before the rest of the data in your frame that is transmitted. The preamble would contain known symbols, so it could be used for multiple purposes, including channel estimation and timing recovery (also known as synchronization).

Once you can recover timing, e.g., from using the received preamble, you can attempt to demodulate your QAM symbols. And besides the phase information that you recover, you also want to avoid any difference in carrier frequency of the signal coming in to the receiver, and the oscillator in the receiver. That's where carrier frequency offset estimation comes into the picture too.


Discrete Fourier Transform and Opposite Convolution Theorem

I'm not sure I understand your argument, but here's a proof which should make the relationship clear:

$$IDFT\left\{\frac{1}{N}(X*Y_N)_k\right\}=\frac{1}{N^2}\sum_{k=0}^{N-1}(X*Y_N)_ke^{i2\pi nk/N}=\\
=\frac{1}{N^2}\sum_{k=0}^{N-1}\sum_{l=0}^{N-1}X_lY_{(k-l)_N}e^{i2\pi nk/N}=
\frac{1}{N^2}\sum_{l=0}^{N-1}X_l\sum_{k=0}^{N-1}Y_{(k-l)_N}e^{i2\pi nk/N}=\\
\frac{1}{N^2}\sum_{l=0}^{N-1}X_le^{i2\pi nl/N}\sum_{k=0}^{N-1}Y_{(k-l)_N}e^{i2\pi n(k-l)/N}=\\=\frac{1}{N}\sum_{l=0}^{N-1}X_le^{i2\pi nl/N}\cdot\frac{1}{N}\sum_{k=0}^{N-1}Y_{k}e^{i2\pi nk/N}=x_ny_n$$

where $(k-l)_N$ means $(k-l)$ modulo $N$. I hope this helps.


wavelet denoising routine for environmental data series

First, a comment - before you denoise, you are basically going to be converting your data from the (time)-domain into the wavelet domain. This is nothing but a series of projections of your data unto user-picked basis functions. (The wavelets).

When you denoise, you will be zeroing out, or shrinking, coefficients which (ostensibly) belong to the noise.

Now, imagine that you picked a wavelet type, and when you transformed your noiseless signal, the very first coefficient was some number, and the rest were zero. This means that you know, that your pure signal, will always give you a non-zero as the first coefficient in the wavelet domain, with other coefficients always being zero.

Now imagine that you transformed your noisy signal. You would probably still see a high value for that first coefficient, (belonging to your true data), but all the other coefficients that you were expecting to be zero have some small value. If you zero them out and inverse transform, you have denoised your signal.

The question you want to answer then becomes, "Which wavelet type, when projected against my data, gives me the sparsest representation in the wavelet domain?" In other words, which wavelet, when used against your (pure, noiseless) data, gives you coefficients that are mostly zero?

Mathematically, what you want is a wavelet that has as many vanishing moments as your data template.

Example: For example, let us say that know that your data template (in your case, the temperature without internal wave motion corruption) is a polynomial of 0th order, (a constant), then you want to use the Haar wavelet, since it has 1 vanishing moment. This means that if you transform your pure-signal into the Haar domain, most of the coefficients will be zero. However when you transform your data-with-internal-wave-corruption, the coefficients that were normally zero now have some value, that you can then clobber and inverse transform.

The same thing can be said for other templates. Let us say that your pure-data was a polynomial with maximum degree of 7. That means then that you want a wavelet that makes my polynomials of maximum degree 7 disappear, or in other words, has 8 vanishing moments. This would make your signal look quite sparse in the transform domain, and thus any other noise of higher order would make it less sparse, and you would be able to denoise by removing those suspicious coefficients that were supposed to be zero.


1) It depends on what you call the standard Kalman filter -- I will call the equations in the picture below to be the "standard Kalman filter". You can easily derive an expression for the Kalman filter where the covariance matrices of the noise processes are time varying in terms of the covariances of the noise processes (and cross-covariance of the noise processes).

Moon and Stirling's Mathematical Methods and Algorithms for Signal Processing derives the Kalman filter under $w(t)$ has mean $0$, covariance $E[w(t) w(s)^T] = Q_t \delta_{ts}$, $v(t)$ is also mean zero with covariance $E[v(t) v(s)^T] = R_t \delta_{ts}$ and $E[w(t) v(s)^T] = M_t \delta_{ts}$ (if you assume $M(t) = 0$, the equations simplify considerably).

For the system model
$ x_{t+1} = A_t x_t + w_t , y_t = C_t x_t + v_t$  [$w_t = w(t), v_t=v(t)$] with $M_t = 0$, the following are the resultant Kalman equations (from Moon & Stirling) to get the predicted state estimate (13.45) and filtered state estimate (13.47):

These equations can be shown to be optimal using either a Bayesian approach (MAP with Gaussianity assumptions) or Linear MMSE approach for the given system model.

2) Convergence of the Kalman filter is a bit of an annoying subject. You can read these notes for more details or check the book Stochastic Systems: Estimation, Identification and Adaptive Control by Pravin Variya and P.R. Kumar or one of the other texts on stochastic systems. Essentially, if the systems are controllable and reachable and yada yada, you can get some statements on convergence properties.


Available Online Recordings for Acoustic Echo Identifiablity

I don't know of any online recordings, but you could try using the Silophone player to create some.  Play your original file, then record what the silophone plays back.
Silophone homepage

The Silophone uses an old grain silo as an echo chamber.  It also plays the echoed audio locally over a speaker outside the silo, so you might want to avoid using anything obnoxious or private or copyrighted.


Downconversion / Demodulation of RF signal to DC

"If the signal spectrum is asymmetric, the original signal spectrum will be corrupted"

If your mixer is a simple multiplication with a single oscillator, then yes, that will be the case.

As Dilip explained, real (RF) signals are always symmetrical in spectrum to the 0Hz axis. That means your original doesn't only exist around 3kHz, it also exists around -3kHz, only mirrored!

I think your image captures that pretty fine:

Now, if you shift your 3kHz signal (black) to 0Hz, the higher spectral peak ends up on the positive frequency side, and the lower on the negative frequency side.

Problem is that for real-valued signals, there's no physical difference between positive and negative frequencies – and thus, your mirrored signal also gets overlayed with this. So you have the small and the big spectral peak overlay - and since they were originally two different parts of your signal, this overlay is not actually just a "shifted copy" of your original signal, but a "corrupted" version.

The only way this wouldn't happen is if the "upper" half of your signal (3-4 kHz) and the "lower" half (2-3 kHz) were already symmetrical in spectrum to its center frequency – because then, the your mirrored spectrum would be identical to the original, and nothing was lost.

But: simply read on. The fact that your text uses negative frequencies definitely works towards showing you how a complex baseband mixer works. In two or three pages, you'll understand what needs to be done to still be


If this is a typical case of a podcast that sounds like a radio broadcast then you could discriminate between music and speech purely based on the "shape" of the spectrum.

I am not sure if there is special software that could handle this but I will provide a simple example using Audacity further below.

The main idea is that a mixed track coming from a variety of different instruments would be expected to occupy a wider range of the spectrum than just someone talking.

Here is a quick test using Audacity:

Just by visual inspection, a set of "noisy bands" appear in the spectrum view. These are the locations where it is either music played along with someone speaking or just music playing.

My settings for this image (Edit->Preferences->Spectrogram): Window-size 4096, window:Blackman, minFreq:0, maxFreq:8kHz, Gain:20, Range:80, freqGain:0, use grayscale.

To switch to this view, load your podcast into Audacity and then click on its track name and select "Spectrogram".

To achieve the same thing in an automatic way, you would have to create a very simple classifier. The classifier would break down the incoming signal into "chunks" and obtain their spectrum. This is roughly what Audacity does to generate each one of the vertical columns that make up the spectrogram of the sound. The classifier would then "scan" these columns from left to right deciding what each column represents ("Music", "Speech", ..., other?). The decision would be based on the relative amplitude of each harmonic in the spectrum. After this, you could decide to drop or silence those frames that have been characterised as "Music".

But all this means writing some form of code in a programming language (e.g. Python). I can expand this response with the simplest of classifiers that could work for this problem but only if that would be seen as useful (?).

Of course, this technique would work with mixed tracks whose full spectrum (including all their possible harmonics) extends well beyond the spectrum of the human voice. There might indeed be some special cases as mentioned earlier by robert bristow-johnson where the discrimination is much more difficult and more complex techniques would be required but for a generic podcast, just basing the decision on the spectrum might do just fine.

By the way, the podcast I used in this example is this one.


Passing a sampled signal through a filter

Not sure what you mean by "wrong".

The Nyquist criteria simply requires you to have "two samples per Hz of bandwidth". It doesn't have to be $[-f_{max},-f_{max}]$, it can be any frequency range that includes at least $2 \cdot f_{max}$ of bandwidth. However, for real signals you need to figure out what to do with the negative frequencies.

For more info google "bandpass sampling".


How do I obtain the frequencies of each value in a 2D FFT and what do they mean?

It's not over your Nyquist limit, because you're doing sampling at your nyquist rate in two independent dimensions and look at the combined result.

Your square-root formula is easily explained: It's pythagoras, matching your combined frequencies in x and y to a single frequency.


What really means stochastic in field of signal processing

Well, getting a bit linguistic, according to the Oxford dictionary:

  stochastic (adj.): Having a random probability distribution or pattern that may be analysed statistically but may not be predicted precisely.

So the definition would be the first one (I don't know where you might have found the second one as you didn't put any source about it). Nevertheless, regarding your specific question, SGD is indeed stochastic. It is an iterative method, but this doesn't mean it is not stochastic. It is iterative and stochastic at the same time: the two concepts are not related to each other.

To put it clearer: if something is blue, couldn't it be big? Well... why not? That's exactly what you are asking here. As size isn't related to colour (absurd example that I think can help here), an approximation can be iterative and stochastic (such as SGD): being one thing doesn't mean it can't be the another one too.


If $N$ is the DFT length, the summation is always over $N$ points. So under this assumption, the first formula doesn't make any sense. Note that the term $e^{-j2\pi nk/N}$ is $N$-periodic in $k$ and $n$, so summing over one period always gives the same value, regardless of the starting index. With the additional assumption that $x[n]$ is periodically continued, i.e., $x[n]=x[n+N]$, you can compute the given sum over any interval of length $N$. A notation emphasizing this fact is
$$X[k]=\sum_{<N>}x[n]e^{-j2\pi nk/N}\tag{1}$$


FIR Impulse response & transfer function

I agree with Matt's comment, however i suspect that whoever gave you the assignment wants you to write down an equation describing the transfer behaviour in time domain for an arbitrary input.

When $x[n] = \delta[n]$, the impulse response of the given system would be

$$h[n] = -\delta[n] + 2\delta[n-1] - 2\delta[n-2] \ .$$

You can find it in the output $y[n]$ by setting the input $x[n]$ to

$$ x[n] = \delta[n] \triangleq \begin{cases}
1 \qquad \text{if } n = 0 \\
0 \qquad \text{if } n \ne 0 \\
\end{cases}$$

Then the output at the end of the diagram is $h[0]=-1$. The $z^{-1}$ elements act as delays of one sample, so the next outputs are $h[1]=2$ and $h[2]=-2$, respectively.

Now to find out what the equation does to general inputs in time domain, write the same thing down in terms of $x[n]$ instead of $x[n]=\delta[n]$.

For any given $x[n]$, the resulting $y[n]$ is produced by the current input times $-1$ and the prior two inputs multiplied with their factors $2$ and $-2$. You can then write down the general solution as

$$y[n] = -x[n] + 2x[n-1] -2x[n-2] \ .$$

Keep in mind however that the term "time domain transfer function" could be misunderstood by some, which is probably the reason your search did not yield any results.