Why is telephone audio sampled at 8 kHz?

If anyone cares to dig I think they will find that before Bell Telephone started multiplexing voice lines they did a lot of research on the frequency content of the human voice.  They originally used test groups to develop the bel audio unit and the power distribution of the voice as well as the sensitivity of the human ear to various frequencies.  They developed a band pass characteristic that peaked around 2.1 KHz and rolled off below 300 and over 3000 HZ.  That gave a good human sounding voice when done correctly.  All that was analog.
AM radio expanded that to 5 KHz to include music that was acceptable to most folks when we were young and had good ears.  Television flyback transformers were designed to run at ~17.5 KHz because there was a magic number for picture reproduction and most folks could not hear the whine.  Single sideband radio was commercialized in the 1960's and needed very sharp cutoff frequencies.  I used radios with filters at 2.1 and 3.1 KHz. 2.1 had some Donald Duck characteristics.  3.1 sounded fine, again with young ears.  Audio bandpass was increased to 20 KHZ or better with FM because the higher carrier frequencies could handle higher bandwidth for better music reproduction.  Stack up some xylophones or bells or other high pitched instruments and they can get enough harmonic energy into the higher frequencies.  OTOH, as was state, most folks cannot hear it.

The bottom line is that anyone claiming they need 20 KHz bandwidth for voice is not paying attention.  3 KHz will do it, 5 will give you some margin.  If it does not sound right then something other than the bandwidth is the issue.

When digital signaling was being developed folks who know figured out that no matter how weird a waveform looked it could be broken down into a set of sine waves.  The harmonic mixing of those waves produced the typical spikey pattern of voice or music.  Lastly, Nyquist did research on the digital sampling rate needed to reproduce a sine wave at a given frequency.  It turns out that it takes 2 samples to make a sine wave so the highest frequency that will be reproduced is half the sample rate.  You want 5 KHz of audio then sample at 10 KHz.  Fine for voice.  You want higher fidelity music than most folks can hear then sample at 40 KHz or so to get 20+ KHz.

One more tidbit is sampling vs bitrate.  If you sample at a given frequency then multiply that by the word length you will get the minimum bitrate needed to produce the desired signal.  Reduce bitrate and the size of the sample word will get cut to meet the new bit rate at a given sample rate.  That is all "lossless" encoding. This is all from memory and trying to find current data.  It's there if somebody looks for cites.  I'm not going to bother as I am getting too old to care.  I just got tired of wading through a lot of obvious mythical issues when I got interested in doing some audio capture.


Hope this is right...
$$h[n]=\delta[n]-a\delta[n-1]$$
$$H(e^{jw})= 1- ae^{-jw} =(1-a\cos(\omega))+ja\sin(\omega)$$
$$|H(e^{j\omega})|= \sqrt{1 + a^2 -2a\cos(\omega)}$$
$$DB_\omega \triangleq 20\log(|H(e^{j\omega})|)$$
$$20\log((1 + a^2 -2a\cos(\omega))^{1/2})= DB_\omega$$
$$1 + a^2 -2a\cos(\omega)= 10^{DB_\omega/10}$$
$$a^2 -2a\cos(\omega) +1-10^{DB_\omega/10}=0$$

with $\omega = 2 \pi \frac{2100}{10000}, \ DB_\omega=1.343 \ \Rightarrow \ a \in \{0.9, -0.4026\}$

$a = -0.4026$ results in a low-pass pre-emphasis whereas $a = 0.9$ results in a high-pass pre-emphasis.


How to sync audio and video signals in MATLAB?

Thanks to helpful comments from @MBaz, I managed to come up with a solution:

we can associate multiple audio samples with a single frame using the VideoFileWriter object. This fact and use-case is missing in the documentation.

First, some stats about the audio and video files. The stereo audio samples are in a 2xN array signal. The video frames are in a cell-array of RGB images, frames.

%% Write audio and video to file
%   Write both audio and video samples into a single video file.
%   Multiple audio samples are matched with one video frame.
%   MATLAB problem: Compression is not possible when audio is included.

% Display A/V stats
fprintf('\nVIDEO STATS:');
fprintf('\nNum. of frames: %d', length(frames));
fprintf('\nFrame rate: %.02f fps', annot.video.frame_rate);
fprintf('\n\nAUDIO STATS:');
fprintf('\nNum. of samples: %d', size(signal,1));
fprintf('\nSampling rate: %.02f Hz\n', infoAudio.SampleRate);

Here's the portion that matches multiple audio samples to a single video frame.
Note that there are 697 video frames and 1336321 audio samples. Some basic fiddling to ensure proper matching. The audio subsampling effect is imperceptible in the output AVI files that I tried.

numAudio = size(signal,1);
numRep = floor(numAudio/length(frames));
numDiff = numAudio - numRep*length(frames); % mismatch

if numDiff
    % if length(frames) does not evenly divide nAudioSamples, then
    % subsample audio to match numRep*length(frames)
    selector = round(linspace(1, numAudio, numRep*length(frames)));
    subSignal = signal(selector, :);
end
assert(numRep*length(frames) == size(subSignal,1));

Finally, we use the VideoFileWriter object available in the Computer Vision System Toolbox of MATLAB to write audio and video to file.

videoFWriter = vision.VideoFileWriter(fullfile(shotPath, 'avclip', ...
                                        [num2str(shotNum) '_av_clip.avi']), ...
                                      'AudioInputPort', true, ...
                                      'FrameRate',  annot.video.frame_rate);
for i = 1:length(frames)
   fprintf('Frame: %d/%d\n', i, length(frames));
   step(videoFWriter, frames{i}, subSignal(numRep*(i-1)+1:numRep*i,:));
end
release(videoFWriter);


The Kalman filter is the least-squares optimal solution for the states of a set of linear dynamic equations under a Gaussian additive noise assumption. Hence, if you have an unforced linear first-order differential equation you can you use the Kalman filter to estimate the state of this differential equation. I'm not sure it will be very useful though, since the convergence rate of the observer has to be much faster than the dynamics of the first-order system in order to capture most of the transient response, which is the only response in an unforced system,  and to converge if the system is unstable (unless you know the initial state, which is rare). High gain amplifies noise, so you might end up with a poor estimate. For such a simple system you might want to use bandwidth-limited numerical differentiation instead to find the inverse of the filter/system to find the state.


Inverse Fourier Transform Dirac impulse with scaled argument

This is just a matter of variable substitution in the integral. The inverse Fourier transform of $\delta(\Omega T)$ is

$$\mathcal{F}^{-1}\big\{\delta(\Omega T)\big\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\Omega T)\,e^{j\Omega t}\,d\Omega\tag{1}$$

Substituting $\alpha=\Omega T$ gives $d\alpha=d\Omega\cdot T$, and $(1)$ becomes

$$\mathcal{F}^{-1}\big\{\delta(\Omega T)\big\}=\frac{1}{2\pi T}\int_{-\infty}^{\infty}\delta(\alpha)\,e^{j\alpha t/T}\,d\alpha=\frac{1}{2\pi T}\tag{2}$$

The last equation is a consequence of the sifting property of the Dirac delta impulse.


I want to perform aec on a wav file and calculate erle

You could use cross correlation for this purpose. For an uncorrelated signal $x(t)$, $\textrm{corr}(x(t), x(t + \Delta t)) = c \delta(\Delta t)$, i.e. an delayed impulse. This can be expanded to comprise multiple echoed copies
$$ \textrm{corr}\left( \sum_{i=0}^{N} a_i x(t + \Delta t_i), \sum_{j=0}^{N} b_j x(t + \Delta t_j) \right)  =
c \sum_{i=0}^{N} \sum_{j=0}^{N} a_i b_j \delta(\Delta t_i - \Delta t_j) $$
I remember from highschool physics books that if a delay is less than a certain amount, associated with the low frequency human hear response, (say 0.1 second), then it is not echo, it is reverberation. More objectively we could say that the audio signal will have a cross correlation that vanishes for lags larger than $L_x$.
You can have a measurement of the echo energy by the sum of the cross correlation with lag longer than $N_e$ from the maximum correlation. The maximum correlation will also provide a reference for synchronization of the two audio signals.
Putting this together you end up with a code like this
def measure_echo(a1, a2 = None):
    corr = signal.correlate(a1, a2 if a2 is not None else a1)
    lags = signal.correlation_lags(len(sig), len(sig_noise))
    # find the peak correlation (synchronize the signals)
    center = np.argmax(np.abs(corr))
    return np.sum(np.abs(corr**2))/np.sum(np.abs(corr[center-Ne+1:center+Ne]))


Calculating the complex signal's average power

Solve the integral:
$$
P = \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} |x(t)|^2 dt
$$
This is usually unwieldy to solve directly; it can be shown that if a signal is periodic, you only need to integrate over the fundamental period $T_0$:
$$
P = \frac{1}{T_0}\int_{-T_0/2}^{T_0/2} |x(t)|^2 dt
$$

EDIT:

Another "trick" is to note that sinusoidal signals of different frequencies are incoherent, so when integrating over long periods of time, it is only necessary to find the power of the individual components because the cross terms go to zero.

This means that
$$
P = P_1 + P_2 + P_3
$$
assuming $x_1$, $x_2$, and $x_3$ are incoherent.  Find the individual powers and sum them.


Understanding Ergodicity and Ensemble Averaging

I will try to explain this practically using MATLAB notation.
Yet before that I must say the ergodic property sometime is limited to a level of moment, namely ergodic in the 1st , 2nd, 3rd moment, etc...
If the process is IID it is promised to be ergodic.

Now, assume we have a function myRandomProcess in MATLAB.
It returns a row vector of length n with the random properties of our random process.

Now, assume we build a matrix A which is $ m $ by $ n $ matrix which each row is a different vector created using the above function.

Now, using the mean and var functions in MATLAB we can compute the properties of the process.

If we apply mean a long the columns (Namely, on each row) we will get $ m $ results.
If those are close enough to the $ n $ results yielded by applying the mean function along the rows (On each column by itself) we can assume the process is ergodic in its first moment.

Now, if $ m $ and $ n $ approaches infinity "close enough" can be small as we want.

Most properties of random process related to its ergodicity can be inferred by looking on its auto correlation function.

Update

Simple MATLAB code to convey the idea:

% Ergodic Random Process

% The Prbability dimension is a long rows (Vertical) and the time dimension
% is a long the columns (Horizontal).

numSamples  = 1000;
numTrials   = 1000;

% Ergodic Process - IID Random Samples

mRandomProcess = randn(numTrials, numSamples);

vTimeAverage = mean(mRandomProcess, 2);
vSampleAverage = mean(mRandomProcess, 1);

vTimeVariance = var(mRandomProcess, 1, 2);
vSampleVariance = var(mRandomProcess, 1, 1);

disp(['Ergodic Process Statistics']);
disp(['1st Moment Ergodicity - ', num2str(mean([vTimeAverage(:) - vSampleAverage(:)] .^ 2))]);
disp(['2nd Moment Ergodicity - ', num2str(mean([vTimeVariance(:) - vSampleVariance(:)] .^ 2))]);

% Un Ergodic Process

% The DC Level is drawed once per random process (Realization)
vDcLevel = (10 * rand(numTrials, 1)) - 5;

mRandomProcess = randn(numTrials, numSamples) + repmat(vDcLevel, [1, numSamples]);

vTimeAverage = mean(mRandomProcess, 2);
vSampleAverage = mean(mRandomProcess, 1);

vTimeVariance = var(mRandomProcess, 1, 2);
vSampleVariance = var(mRandomProcess, 1, 1);

disp(['Un Ergodic Process Statistics']);
disp(['1st Moment Ergodicity - ', num2str(mean([vTimeAverage(:) - vSampleAverage(:)] .^ 2))]);
disp(['2nd Moment Ergodicity - ', num2str(mean([vTimeVariance(:) - vSampleVariance(:)] .^ 2))]);

As you'll be able to see, the first process is Ergodic (As any IID process), as it is equivalent to calculate its moment along the time dimension or the sample dimension.

On the other hand, the second one isn't Ergodic.
Since its mean along the rows is zero (Drawn from a uniformly random variable [-5, 5]) yet along its columns its mean is the value of the specific realization of the uniform variable.

I hope this clarifies it.


Impulse Response to Frequency Response in Octave

The Frequency Response of a system is the Fourier Transform of the Impulse Response of the system.

Since we're in the real world and we have finite number of samples observed over finite time interval we use the Discrete Fourier Transform (DFT).
In Ocatve and MATLAB the DFT is implemented in efficient way using the Fast Fourier Transform (FFT).

In order to calculate it on MATLAB / Octave you can do the following:

% Assuming Impulse Response is given by 'vImpulseResponse'

vFreqResponse = fft(vImpulseResponse);

figure();
plot(abs(vFreqResponse));

Enjoy.


Are there accepted ranges of "high", "medium" and "low" visual spatial frequency?

The short answer is yes(1), no(2).

The quantification of the power spectrum is performed via the Fourier Transform over an image. An image is a matrix that describes the distribution of reflectance (or radiation more generally) across the field of view of the camera. It therefore has only a remote "relationship" with what is actually depicted in the image (the physical scene). Image resolution and spatial resolution are two different things. You can use a 50MP camera at the end of a telescope or at the end of a microscope. Both images will have the same image resolution but they will be "looking at" different spatial resolutions (actual space being sampled).

In other words, small objects correspond to high frequencies and large objects correspond to low frequencies but what is large and small depends on the context of the image (and the characteristics of the imaging equipment to an extent). Consequently, what is "high", "medium" and "low" spatial frequency can only be expressed in relative terms and either with a lot of assumptions or constraints.

The approach of dividing the available bandwidth into three bands and estimating the average power in each is, of course, reasonable.

An improvement would be to take a set of representative images, average their spectra and look at the actual bandwidth they occupy and its shape with a reference to the actual objects they represent to then decide at a different possible assignment between the low, med, high characterisations. But this implies specific scenes with a specific field of view, orientation, resolution and any other factor that affects the recorded spatial frequencies.

In case it is of any additional help, there is a concept in image histograms called "Adaptive Binning" where, the width of a given histogram bin is not fixed. The width is decided after an optimisation process over a criterion and it is image specific. For example, one criterion could be that each bin should represent at least N pixels. Although originally not defined over spectra, this adaptive concept might be useful here as well but applied over FFT bins. For more information please see this link.

Hope this helps.


For linear IIR digital filter, what happens for negative frequencies?

It may be easier to start with analog signals. If the filter's impulse response $h(t)$ is real, then its frequency response $H(f)$ is conjugate symmetric: $H(-f)=H^*(f)$, where $(\cdot)^*$ indicates complex conjugate. Since $H(f)=|H(f)|e^{j\angle H(f)}$, then conjugate symmetry implies that $|H(f)|=|H(-f)|$ and $\angle H(f) = - \angle H(-f)$. This is true for any kind of filter, FIR or IIR.

So, if the filter input $x(t)$ is $$x(t)=A_0\cos(2\pi f_0t+\phi_0)=\frac{A_0}{2}\left(e^{j2\pi f_0t}e^{j\phi_0}+e^{-j2\pi f_0t}e^{-j\phi_0}\right),$$ then the filter output is

$\begin{align*}
y(t) &= \frac{|H(f_0)|A_0}{2}\left(e^{j2\pi f_0t}e^{j(\phi_0+\angle H(f_0))}+e^{-j2\pi f_0t}e^{-j(\phi_0+\angle H(f_0)}\right) \\
&= |H(f_0)|A_0\cos\left(2\pi f_ot+\phi_0+\angle H(f_0)\right).
\end{align*}$

As you can see, the amplitude of negative frequencies are affected exactly the same as positive frequencies. The phases are affected in a similar fashion, but with phase of opposite sign for negative frequencies. Note that this is necessary to preserve the symmetries that real signals have in the frequency domain.

In the case of discrete signals and filters, the conclusion regarding negative frequencies is the same. Finding the actual gain and phase at a particular frequency is slightly more difficult, though. The reason is that the DFT of the filter's impulse response defines the filter's gain and phase only at a finite number of frequencies. In general, the filter's gain and phase at an arbitrary frequency $f_0$ can be found by evaluating the filter's transfer function at $z=e^{j2\pi f_0}\,$.


You're conflating two different definitions of orthogonal.

Orthogonal in the OFDM sense means that the subcarriers are mutually orthogonal to one another from the perspective of a correlator-based receiver (the optimum receiver for the AWGN channel).
An orthogonal matrix $\mathbf{A}$ is one whose columns are all unit vectors that are jointly orthogonal to one another, such that $\mathbf{AA^H} = I$.

These two concepts aren't related; there is no constraint on the modulation applied to the various subcarriers in an OFDM symbol. The orthogonality is imposed via careful choice of the symbol period and subcarrier spacing, typically implemented using DFTs since it conveniently provides the desired time/frequency grid.


Load and Play Wave File (`.wav`) on MATLAB

When you load Audio file to MATLAB you basically load array of numbers.
When you play Audio there is a parameter which is needed - Sampling Rate.
This is what tells the player how fast to play the data.

Think if you load a file with 1000 samples. How long should is the audio?
Unless you supply some information, the player can not know (MATLAB just use a default which doesn't match your data).

If you look at the documentation of the sound() function you'd see there is an option to use it as sound(audioSignal, samplingRate).

In order to get the Sampling Rate of a Wave file use: [audioSignal, samplingRate] = audioread(filePath). Read more at audioread().


Estimate a Transfer Function from ARX Models vs. ARIMAX?

Hi: I don't know about BJ or OE but, as far as the other three go, it REALLY MATTERS which one you choose because ARX has no MA or differencing terms, ARMAX has no differencing terms and ARIMAX has both. But, even within any of the three models, there are major choices to be made as far as A) the number of lags of y, the number of lags of  X  ( and, in the ARDL case, the number of X's). This is an EE forum so  I'm only on it because I'm trying to learn some EE. My background is more econometrics so my suggestion is that you should probably look at a good econometrics textbook that discusses model selection given the question you asked. There is the top down method of hendry, the schwert criterion, AIC etc.  Luktepohl's text probably discusses model selection in these types of models but focuses on AR and ARI and not so much the MA terms. Also, it sounds like the response y_t is a scalar but, in case it isn't, you want to look at VARS which luktepohl covers in great detail. Good luck. Oh,  Also, you may want to look at the koyck distributed lag which is a special ARMAX model.


Shannon interpolation formula for downsampled data with an "almost ideal" low pass filter

I don't get your downsample step when you downsampled by factor $M$.

Let me go from scratch with the spectrum visualization below, with time domain, continuous frequency domain and discrete frequency domain from left to right.

When we reduce the sampling frequency by a factor $k$, the signal spectrum is copied to new replicas at $f_s/k$. The discrete spectrum is a snapshot of continuos spectrum at $[-f_s, f_s]$ and $[-f_{s,down} = f_s/k, f_{s,down} = f_s/k]$, so it is expanded by the factor $k$.

The DFT works on the discrete frequency domain. The downsample $k$ must be chosen so that the expanded spectrum (in discrete frequency, or the new replicas in continuous version) does not overlap with their copies centered at $-2\pi$ and $2\pi$, or $-1/2$ and $1/2$ if you take discrete instantenous frequency by dividing $2\pi$; otherwise aliasing happens and the reconstruction of $x(t)$ is impossible. If there is not alias, a low pass filter to take the spectrum part centered at 0 is suffice. This filter does not need to be "ideal", just to be sure that the filter takes only the center part.

In your calculation, if I understand well, you are talking about the discrete frequency domain and the original spectrum takes $1/M$ the normalized band (you said $f$ is in $[-1/2, 1/2]$ and your signal has support in $[-1/2M, 1/2M]$). In this case your downsample factor $k$ must be less than $M/2$ and the cutoff frequency of your ideal filter must be at $1/4$ if $k=M/2$ (and yes, in this case $k=M/2$, we need "ideal filter" assumption).


When is a discrete time transfer function unrealizable?

You give a transfer function in $z^{-1}$, but for the purposes of making a state-space model, you need a transfer function in $z$ (more on that later).  So your $H$ becomes
$$H(z) = \frac{b_0 z + b_1}{z^2 + a_1 z + a_2}$$
When you invert that it's improper, which shows that you need to know information one time step into the future.  That's not possible (or, if it is, you should take up life as a stock broker).

(Note that you don't need state space to know this: you can just write the corresponding difference equation for your $H^{-1}$ and you'll see that it's asking for future values of the input to affect the current value of the output.)

Appendix

The reasoning for doing it in $z$ is because the time-domain state-space realization is

$$\begin{split} x_k = A x_{k-1} + B u_k \\ y_k = C x_{k-1} + D u_k\end{split}$$

For $D = 0$, the resulting transfer function is strictly proper in $z$ (but not in $z^{-1}$); for $D \ne 0$ the resulting transfer function in $z$ has equal orders in the numerator and denominator.


How sampling aperiodic signal will result in periodic repetitions of the same

sampling guarantees periodicity

Not exactly. Sampling in one domain guarantees periodicity in the other domain. So sampling in time creates a periodic spectrum and sampling in frequency creates a periodic time signal. Or the other one around: a periodic time signal has a discrete spectrum and a periodic spectrum has a discrete time signal.

That's why there are four different "types" of Fourier Transform. One for each permutations of continuous vs discrete in each domain. See for example http://fourier.eng.hmc.edu/e101/lectures/handout4/node3.html or https://www.dspguide.com/ch8/1.htm


How to find ROC of system when input is two sided and output is one-sided

The transfer function is

$$H(z)=\frac{Y(z)}{X(z)}=-\frac18\frac{z-3}{\left(z-\frac14\right)(z-\frac{10}{3})}\tag{1}$$

It has three possible ROCs:

$|z|<\frac14$ (left-sided, unstable)
$\frac14<|z|<\frac{10}{3}$ (two-sided, stable)
$|z|>\frac{10}{3}$ (right-sided, unstable)

Only the second ROC, corresponding to a two-sided and stable impulse response, overlaps with the ROC of $X(z)$. The intersection, and, consequently, the ROC of $Y(z)$ would be $\frac13<|z|<3$ (equal to the ROC of $X(z)$). However, the pole at $z=3$ is removed by pole-zero cancellation, which results in the final ROC $|z|>\frac13$ for $Y(z)$.


How is a PCM byte converted to a voltage level?

There are quite a few different DAC topologies with different properties. The typical tradeoff is between resolution and speed and power usage -- higher speed means lower resolution, and higher resolution means higher power usage.

The simplest DAC is arguably the string DAC or Kelvin divider. If you consider the 3-bit DAC above, you can see that every possible value representable by 3 bits has a corresponding voltage level in the string. As indicated, you decode the 3-bit word to a switch a voltage level to the output, turning off all the others. For example, if you have the word 010, you turn on the 3rd switch, and if you have 101, you turn on the 6th.

You can extend this string to have $2^{16} = 65 536$ levels, but this make a very large circuit, and is not feasible in practice. The more compact R-2R ladder is then a better option, although a pure 16-bit R-2R ladders will have very high requirements for component tolerances. A combined string and R-2R topology is quite common for 16-bit.

For audio, the so-called $\Delta-\Sigma$ modulator is quite popular. In this case the output is switched, and averaging in the form of a low-pass filter is used to generate the desired levels.


Does gray coding bring improvement to SER of QPSK, compared with natural coding?

The central point as crystallized in the comments is that you think SER is derived from BER.

It's not!

SER is short for "symbol error rate", and that's it. That's the underlying thing.

A symbol error is when one symbol was sent, and a different one was received. It doesn't matter whether the symbol you've sent was "labeled" 00, or 11, or 10; or whether the received symbol was defined to be 00, 11 or 10, or yoghurt for that matter; what matters is that not the symbol that the sender sent was decided for in the receiver.

Thus, the SER doesn't care about what binary values are assigned to the symbols. (and that assignment can be Gray, or something else.)

By the way, I'm not sure what you mean with "natural encoding": I know that term from one-dimensional encoding problems (e.g. ADC values), where natural is just binary counting up:

value    |  0  1  2  3
---------+----------------
encoding | 00 01 10 11

If you extend it two a second dimension by just forming coordinate vectors, you get gray encoding for QPSK:

    \ dim0  0   1
      -------------
dim1 |
     |
   0 |     00  01
     |
   1 |     10  11

So, I'd argue Gray and natural are the same here, and using anything different would just be ill-advised .


Phase changes of BPSK modulated carrier after moving to baseband

There is no change - the complex envelope of the signal (magnitude and phase) at any given carrier frequency is no different from when that carrier frequency = 0 (DC, at baseband). To hold this equivalence the baseband signal must be complex (often represented with I and Q data paths) unless the signal is completely on the real axis... for example If the BPSK Signal at the carrier has a 10 degree phase offset (such that the carrier frequency is toggling between 10 degrees and 190 degrees), the baseband signal will do the same as long as it is allowed to be complex.

For this reason it usually doesn’t make sense to actually model a sinusoidal carrier when doing simulations but work directly with the complex baseband signal.

The further details given on upconversion and downconversion in this post may help provide further insight:
Frequency shifting of a quadrature mixed signal


How to determine the filter type of a discrete frequency spectrum

The $k_1$ indexes you get from the FFT correspond to frequencies 0 to $2\pi$ (actually, 1 sample short of $2\pi$). Note that frequencies in the range $\pi$ to $2\pi$ also correspond to the negative half of the specturm ($-\pi$ to 0).

You can get a nice visualisation with frequency 0 in the middle by applying fftshift to the result (Octave or Matlab assumed).

You can also check freqz, which plots a nice frequency and phase repsonse (try freqz([1 2 1], [1])). Or try plot(abs(freqz([1 2 1],[1]))) to get a linear plot rather than log.

What you are getting is a low-pass filter (gain 4 at low frequency, 0 at high frequency).

In VERY general terms: You can know beforehand it is a LPF because the impulse response is a single bump. A high-pass filter will have some sort of antisymetric shape (like [1 2 -2 -1]), and a band-pass or band-reject will have some oscillations in the impulse response.


how to detect corner in a digital image using MATLAB?

In the case of video, the best approach to use is the ShiTomasi Corner Detector. (detectMinEigenFeatures). Now, if you want to detect corners that contains some invariant featrues, you should use SURF, BRISK, ORG, AKAZE, FREAK, etc...

Example:

im = rgb2gray(imread('image.jpg'));
[r c] = corner_ST(im,20);


Answers in textbooks are not necessarily right.. I have got the same answer as you ;)

If you are not sure that your method is right, you can always test digital systems on your PC (matlab, scilab, python, octave, or almost in any language )


The OP is trying to do what we would refer to as channel estimation. I am linking two posts below where this is all worked out in greater detail but do want to make the point first that the choice of "sounding pattern" is important. The channel can only be estimated where there is signal energy in the frequency domain, so the sine sweep is an excellent choice as long as the sweep is over all the frequencies of interest. Other sounding patterns commonly used are pseudo-random noise generated with PN sequences (from linear feedback shift registers: LFSR).
An optimized FFT chirp for channel estimation:
How can I plot the frequency response on a bode diagram with Fast Fourier Transform?
Using Wiener-Hopf equations for least-squares channel estimation:
Compensating Loudspeaker frequency response in an audio signal
Related posts on PN sequences as an alternate sounding pattern:
What are the advantages and disadvantages of PN sequence over Walsh code
Sequence Length for a Linear Feedback Shift Register


1) Do I fail to understand the question? The allowed frequency range of a bandpass filter goes from Fpass1 to Fpass2. The range between Fstop1 and Fpass1 and from Fpass2 to Fstop2 is the transition range where the filter goes from fully attenuated to fully "open" so there shouldn't be any strong interferers in this range.

2) Your band select filter before the down-mixer. It usually has your carrier frequency in the middle of Fpass1 and Fpass2 and is as wide as your expected signal is.

3) 999000


What Is the Derivative of the Function $ s \left( t \right) = \left( 1 - {e}^{\frac{-t}{RC}} \right) u \left( t \right) $?

By the Product Rule (With abuse of the derivative of Unit Step):

$$\begin{align*}
\frac{d}{dt} \left[ \left( 1 - {e}^{\frac{-t}{R C}} \right) u \left( t \right) \right] & = \frac{d}{dt} \left[ \left( 1 - {e}^{\frac{-t}{R C}} \right) \right] u \left( t \right) + \left( 1 - {e}^{\frac{-t}{R C}} \right) \frac{d}{dt} u \left( t \right) \\
& = \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) + \left( 1 - {e}^{\frac{-t}{R C}} \right) \delta \left( t \right)
\end{align*}$$

By the properties of Dirac's Delta $ f \left( t \right) \delta \left( t \right) = f \left( 0 \right) \delta \left( t \right) $ we get:

$$ \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) + \left( 1 - {e}^{\frac{-t}{R C}} \right) \delta \left( t \right) = \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) + 0 \delta \left( t \right) = \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) $$

Remark
The above is with some abuse of Math (Derivative of Unit Step as Delta Function, Zero multiplied by Delta is zero, etc...) which is used (Again, it is abuse but usually gets to the right place).
If I remember correctly, even this intuitive abuse is used in Desor.


Convolution of $\sin(\omega t)$ and $\cos(\omega t)$?

This is a theoretical question without much practical interest but still it can be nice to check the results and investigate if intuition still holds.

First of all, if the convolution is regarded an LTI operation between an input $x(t)=\cos(\omega_0 t)$ and a system $y(t)=h(t)=\sin(\omega_o t)$ then it's immediately obvious that since $h(t)$ is an unstable system, the output can be unbounded even for a bounded input signal.

Furthermore assuming that Fourier transforms of the input, the system and the output exists (the integrals converge!) then we can assert the property that: $$ x(t) \star y(t) \longleftrightarrow X(\omega)Y(\omega) $$

When the sinusodal signals $x(t)=\cos(\omega_0 t)$ and $y(t)=\sin(\omega_0 t)$ are considered, it's obvious that their Fourier integrals do not converge. Hence their formal Fourier transforms do not exist. The solution is the acceptance of the generalised impulse function to represent the Fourier transform of the sinusodial signals; i.e., $$ x(t)=\cos(\omega_0 t) \longleftrightarrow X(\omega) = \pi [\delta(\omega-\omega_0) + \delta(\omega+\omega_0)]$$ and $$ y(t)=\sin(\omega_0 t) \longleftrightarrow Y(\omega) = \frac{\pi}{j} [\delta(\omega-\omega_0) - \delta(\omega+\omega_0)]$$

Then we consider that the above theorem still holds for the impulse functions as well and apply it here:
$$ z(t)=\sin(\omega_0 t) \star \cos(\omega_0 t) \longleftrightarrow Z(\omega)=\frac{\pi}{j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)] \cdot \pi[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)] $$

Using the impulse sifting property; $$f(x)\delta(x-a) = f(a)\delta(x-a)$$
we shall perform th multiplications as follows:
$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(\omega-\omega_0) \delta(\omega-\omega_0) + \delta(\omega-\omega_0) \delta(\omega+\omega_0) - \delta(\omega+\omega_0) \delta(\omega-\omega_0) - \delta(\omega+\omega_0) \delta(\omega+\omega_0)] $$

$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(0) \delta(\omega-\omega_0) + \delta(-2\omega_0) \delta(\omega+\omega_0) - \delta(2\omega_0) \delta(\omega-\omega_0) - \delta(0) \delta(\omega-\omega_0)] $$

Now noting that $\delta(2\omega_0)=\delta(-2\omega_0)= 0$ those terms cancel and $\delta(0)=\infty$ terms remain, hence

$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(0) \delta(\omega-\omega_0) - \delta(0) \delta(\omega-\omega_0)] $$

$$ Z(\omega)= \pi \delta(0) \frac{\pi}{j} [ \delta(\omega-\omega_0) - \delta(\omega-\omega_0)] $$

Which is recognized as the Fourier transform of the infinite amplitude sine wave as: $$ \boxed{ z(t) = \pi \delta(0) \sin(\omega_0 t) } $$

The conclusion is that the convolution between $x(t)=\sin(\omega_0 t)$ and $y(t)=\sin(\omega_0 t)$ produces an infinite amplitude sinudoidal wave of the same frequency $\omega_0$.

The time domain verification is as follows:

$$x(t) \star y(t) = \int_{-\infty}^{\infty} x(t-\tau)y(\tau) d\tau \leftrightarrow z(t) = \int_{-\infty}^{\infty} \cos(\omega_0(t-\tau)) \sin(\omega_0 \tau) d\tau$$

Using the trigonometric identity of $$\cos(x)\sin(y) = 0.5[\sin(y+x) + \sin(y-x)]$$ we can break the integral into two, where $x = \omega_0(t-\tau)$ and $y=\omega_0\tau$ , hence
$$z(t) = 0.5 \int_{-\infty}^{\infty} \sin(\omega_0(t-\tau)+\omega_0 \tau) + \sin(\omega_0 \tau - \omega_0(t-\tau) ) d\tau$$

$$z(t) = 0.5 \int_{-\infty}^{\infty} \sin(\omega_0 t)d\tau + 0.5 \int_{-\infty}^{\infty} \sin(2\omega_0 \tau - \omega_0 t) d\tau$$

Now the first integral becomes $$ 0.5 \sin(\omega_0 t) \int_{-\infty}^{\infty} d\tau$$ while the second integral can be shown to be zero after a suitable change of variables, assuming for a given (fixed) t, $\phi = 2\omega_0 \tau - \omega_0 t$ and $d\phi = 2\omega_0 d\tau$ then the second integral becomes $ \frac{1}{2\omega_0} \int_{-\infty}^{\infty} \sin(\phi) d\phi = 0$.

Theferore the result of the convolution is:
$$z(t) = \left( 0.5 \int_{-\infty}^{\infty} 1 d\tau \right)\sin(\omega_0 t) $$
Note that the integral that weights the sine wave has infinte value, moreover it can be shown that by the forward and inverse Fourier transform pairs one can recognize the integral as the forward Fouier transform $H(\omega)$ (evaluated at $w=0$, $H(0)$) of the signal $x(t)=1$ which is:

$$\mathcal{F} \{ 1\} = \int_{-\infty}^{\infty} 1 e^{-j \omega t} dt  \equiv 2\pi \delta(\omega) $$ and therefore setting $w=0$ yields $$\int_{-\infty}^{\infty} 1 dt  \equiv 2\pi \delta(0) $$

Finally plugging this into the result yields; $$\boxed{ z(t) = \pi \delta(0) \sin(\omega_0 t) }$$ which is the same as the result obtained from the Fourier method.


the problems seems to be with the white border of the image you are using. remove it manually and test again.


Sampling frequency and oversampling factor

The rule you must remember when deciding on your sampling frequency is: it must be at least double the highest frequency signal that will ever appear in your simulation.

Consider your baseband signal. The highest frequency signal is the half-sine pulse (its bandwidth is infinite in theory). We know, from the Fourier series of a full-wave rectified sine, that its spectrum decays with frequency as $$\frac{1}{\pi(4n^2-1)}.$$ I would go with $n=3$, which corresponds to a frequency of $3f_0$. If your pulse rate is 2 Mpulses per second, then $f_0=2$ MHz, the pulse bandwidth is 6 MHz, and your sampling rate must be 12 MHz.

Regarding passband: If you must simulate this system at a high carrier frequency, you know that the bandwidth now will be 2.4 GHz plus 6 MHz, or 2.406 GHz, so your sampling frequency needs to be double that. However, before doing this, please learn about the complex envelope. There is nothing to be gained by simulating the passband system, and you can lose a lot of time.


What is the correct solution for Fourier transform of unit step signal?

As Matt said, the second and third definition are the same except for the part with impulse. The impulse ($\pi\delta(\omega)$) account for the DC value of $u[n]$. Without that term (i.e., the second definition) it is actually the FT of $v[n]=\frac{1}{2}\operatorname{sgn}[n]$. We have $u[n]=v[n]+\frac{1}{2}$. And hence the FT of $u[n]$ has the additional term to account for the addition of $\frac{1}{2}$. Also, the discrete time FT (or DTFT) of $u[n]$ is properly written as $U(e^{j\omega})$.

The first definition, $U(j\omega)$ is the "continuous time" FT (or CTFT) of $u(t)$ (not $u[n]$) and hence different from the other two definitions.


The frequency-domain spectrum obtained using MATLAB's FFT function seems inaccurate

The FFT produces values for both negative and positive frequencies.

Only plot the first half of the  FFT output.


What's the number of data points in a Welch-based PSD?

And that's absolutely correct. As you mentioned, procedure can be described as follows:

Take long signal of length $N$ and slice it with some overlap and window into segments of length $L$.
For each of these segments calculate the DFT. Assuming no zero padding, you will get spectrum with $L$ points - same as segment length.
Because signals are real, then we just care about first half of spectra (squred magnitude to be correct), thus take only $L/2$ points.
Do this procedure for each segment and you will end up with lot's of vectors (dependent of overlapping and signal length $N$) of length $L/2$.
Now you need to average all the spectra. This is simply average for each frequency bin across number of segments and PSD is calculated.

Let's assume you will end up with 50 segments of length $L$, then you must average 50 spectra of length $L/2$. The longer is your signal, more segments will be averaged (better estimate), but obviously their length is not changing - you've chosen it in the beginning.


Detecting Pattern from Signal Data by Gaussian Mixture Model?

If the data is cyclic by its nature the best thing would work using its spectrum.

You can easily build a system which checks sub set of data to verify periodic and the once you establish your groups checking the affinity of new data is easy - add it to each series. It should belong to the one creates less spread in the frequency (Namely, it follows the cycles).


How to get coordinates(x(time),y(frequency)) of spectral peaks from spectrogram?

Under the hood, your spectrogram consists of multiple, consecutive FFT frames (of length N) and each FFT frame contains N/2 frequency bins (even if you have N frequency bins, the upper part is discarded).

An FFT frame denotes a temporal position from the beginning of a signal and this temporal position can be computed as

time_pos = frame_number *  N / sampling_rate; //frame_number runs from 0......N-1

As an example, if your signal is sampled at 8000 Hz and N is 1024, the 10th frame will be at 10 * 1024 / 8000  = 1.28 sec from the beginning of the signal.

Just like every FFT frame has a fixed length (N samples long), each frequency bin has a fixed width in Hz (frequency resolution). This frequency resolution is given by the sampling rate divided by the frame/fft length so

freq_res = sampling_rate / N;

In order to find the frequency for a particular bin, simply multiply the bin index by the frequency resolution.

freq_10bin = 10 * freq_res;

In order to get the coordinates of the spectral peaks, you should iterate through all the frames and find those bins that have a value above some threshold thus indicating a peak.
The frame indexes and bin indexes respectively are then used to compute the temporal and frequency coordinates for every peak (and the corresponding displacements).


Conceptual Question from Signal Processing - Impulse Response and AR Coefficients

Impulse Response is basically the FIR coefficients of the system.
Namely, a system $ H $ with an impulse response given by $ f [n] $ and a Filter $ F $ with an FIR representation of $ {f[0], f[1], \cdots, f[n]} $ are equivalent.

Now, systems with Feedback are equivalent of both FIR and IIR (AR) filters.
But given infinite length of FIR model any LTI system can be represented by FIR coefficients.

The relation could be easily displayed by the Laplace Transform of the system.


Video processing to detect a logo in frames

If you know where the logo is then you can detect its presence, fairly quickly, with cross-correlation...with a little bit more work.

Detecting the presence of "things" using cross-correlation is a way to do template matching, where we are trying to find a specific pattern in a longer "duration" signal.

In this case, the logo is the pattern and the longer "duration" (or surface in the case of an image) is the whole frame.

The question now becomes, how do you recover the pattern? How do you recover the logo, so that you can then search for that particular target?

There are two ways to do this:

It is already available: In which case, you just download a copy of the logo, resize it to the approximate size you expect the logo to be in the frame and then use cross-correlation to find it.
It is not available and it has to be recovered by the frame itself: To recover the logo, use averaging. Take a large number of frames (say 30-50 seconds of video), and average their pixel values. The key idea here is that the logo does not "move" from its position. Therefore, the pixels that correspond to the logo will have a stable value while the pixels that correspond to the background will be taking "random" values from whatever the camera happens to pan to. When you average, the stable values will retain their magnitude and the random values will have considerable lower values. The reason I am telling you to average such a long duration is because the second image is showing a frame from a talk-show. In that one, the camera, the subject or its background are not expected to move too much, therefore, it is likely that if you average for a short duration you might "pick up" a bit of the background too. In that case, you might have to recover the logo from a different part of the video. The bottom line is that this process will give you the form of the logo. You might have to apply thresholding to get rid of those "background" pixels to get a proper "mask" image.

Once you have that small "logo" image, you can then pick up a random frame from your video feed, crop it around the approximate location of the logo and then take that part of the frame and cross-correlate it with the logo. If the logo IS present in the frame, you will get a much stronger cross-correlation peak at its approximate centre than if the logo is not present. You can then apply a threshold to detect this strong peak. If the peak is above that threshold then the logo is present, if not, then the logo is not present. That's the concept and different platforms will have different implementations (e.g. here is the same thing in Python).

Hope this helps.


Why are recursives methods useful for FIR filter design?

Sure, FIR filters are "simple" linear objects. But a typical FIR filter design task looks like this:

  Find a set of numbers (the filter coefficients) such that the
  magnitude of their Discrete Fourier Transform is as close as possible to a
  target response. "As close as possible" being defined as "by
  minimizing the $L_\infty$ norm of the difference between ideal and actual
  response over some frequency ranges".

This is a non-linear, sometimes non-convex optimization problem and there is no universal, trivial procedure to solve those.

We have to use procedures like the Parks-McClellan algorithm because the real-world engineering constraints we want to impose when designing FIR filters do not translate into simple linear or quadratic constraints. For example, a least-square error criterion (L2 norm often leads to more tractable procedures than other norms...) would not be an acceptable choice for engineering applications, because it could allow solutions that are close to the required response on average, but with an outlier value.


Why does block LMS have the same performance as LMS?

here is a short and little (a.k.a. "consise but terse") derivation of the LMS and normalized LMS adaptive filter.

once the LMS has converged on a reasonably stable equilibrium for $h_n[k]$, they won't move around that much.  then it doesn't matter so much how long the block is.  and the only difference between block-LMS and the plain-old ordinary LMS is the block size.  (the block size for the latter is 1.)


Audio transfer function modelling with pink noise + signal

your generated pink noise is just a signal.  it's not noise because you know what it is.

you can toss whatever broadbanded signal you want into your linear, time-invariant system.  call that "$x[n]$".
you can measure the output of that LTI system, call that "$y[n]$", given the known input.
you can FFT both the input and the output.  call those FFT results "$X[k]$" and "$Y[k]$", respectively.
if the input was broadbanded, $X[k] \ne 0$ for all $k$.  so you can divide $Y[k]$ with $X[k]$.
the result: $$ H[k]=\frac{Y[k]}{X[k]} $$ is your system response.  this is called the "Dual-channel FFT" method of measuring system response. if there is a delay, you will see a linear trend to the unwrapped angle of $H[k]$. the negative of the slope of that linear trend is proportional to the delay.

you can also use cross-correlation between the input and output to discover the prominent delay of the system.  the definition of that is

$$ R_{xy}[m] = \sum_n x[n] y[n+m] $$

i'm deliberately playing fast-and-loose with the limits of the summation.  the value of $m$ that maximizes $R_{xy}[m]$ is a good candidate for the delay of the system.


What's the best way to compare the same object from different "crops" or angles within multiple images?

There are many ways to achieve this. You could use an image based, direct similarity metric such as Sum of Absolute Distances, Sum of Squared Distances or Normalized Cross Correlation. For a general overview look here. These approaches generalize to template matching where a convolutional procedure is conducted, such as the one here.

Template matching methods could tolerate rigid transformations (rotation, translation) and scaling. However, for the case of affine deformations (e.g. the second image set in the question), they would not suffice. It is then better to benefit from feature based matching techniques, where some repeatable keypoints are first extracted from the image, and invariant descriptors of the local neighborhoods are computed and store. This is done for both images. Finally the descriptor are compared using a nearest neighbor search. SIFT is one of those robust descriptors.

Note that, while in the first scenario, the template transformation can be defined simply by a translation, the second case is not like that. You should either solve for affine parameters, or compute a general perspective homography to represent the pose of the second image w.r.t the first one.

The problem is very general, and often times analyzed under the roof of object detection. You could review this literature and choose a method suited to your problems. OpenCV implements most of the techniques, I described in this answer.


Pixel Binning: Effect on SNR for Hardware vs. Software Binning

If you have access to the raw pixel data, you can implement binning in software too. They should be identical upto the step-size of the ADC, by which the quantized analog sum of pixels may differ from sum of quantized pixels.

Hardware binning has the advantage that it has to transfer less number of pixels, hence can operate at high fps compared to raw data transfer.


What is the advantage of performing the FSK using IQ modulation?

You only need one RF oscillator operating at a single frequency, instead of having 2 (or more in the case of M-ary FSK) oscillators operating at separate frequencies for each bit/symbol.
Since you only have one oscillator, you don't have to worry about discontinuities in the phase of the transmitted RF signal, like this:

Sharp discontinuities contain high frequency which widens the bandwidth of the signal. With I/Q modulation you'll have a much easier time producing a continuous-phase FSK signal (CPFSK).

  How to demodulate the FSK signal?

First lets assume we are trying to demodulate a simple FSK signal where the center frequency is $f_c$, a '1' is represented by the frequency $f_c+\Delta f$, and a '0' is represented by the frequency $f_c-\Delta f$.

Before IQ sampling, you tune your RF front end receiver to $f_c$. Then when you sample you have a baseband I/Q representation of the signal, centered around 0 Hz. The only thing to do is figure out if you have a positive frequency (+$\Delta f$) or negative frequency (-$\Delta f$). The easiest way to do that is to monitor the change in phase of the baseband IQ signal. A positive phase change corresponds to a positive frequency, and a negative phase change corresponds to a negative frequency.

So monitor the phase (the angle of $I + jQ$) and watch how it changes over one symbol period. Positive change = '1', negative change = '0'.

  Does the FSK signal look triangular wave or sinewave?

If you're talking about the RF signal, it looks like a sinewave of varying frequency.

If you're talking about baseband IQ signal, it also looks sinusoidal, but with a bunch of phase shifts, like this (blue = I, orange = Q):

In this particular example, a '1' corresponds to a phase change of $+\pi/2$, and a '0' corresponds to a phase change of $-\pi/2$.


How can one improve the robustness of adaptive beamformers to signal mismatches?

One approach that I'm aware of is to "load" the diagonal of $\widehat{\mathbf{R}}$ by some value $\sigma_d^2$ in order to stabilize the covariance matrix and make the adaptive vectors more robust to signal mismatches. This is often done in matched field processing and is referred to as diagonal loading or white noise constraint beamforming.

In other words, the robustified adaptive weight vectors are calculated as

$$\widehat{\mathbf{w}}_\theta=\frac{\left(\widehat{\mathbf{R}}+\sigma_d^2\mathbf{I}\right)^{-1} \mathbf{a}_\theta}{\mathbf{a}^H_\theta \left(\widehat{\mathbf{R}}+\sigma_d^2\mathbf{I}\right)^{-1} \mathbf{a}_\theta}$$

and the output beam power is now computed as

$$\widehat{B}_{mvdr}(\theta)=\widehat{\mathbf{w}}_\theta^H \left(\widehat{\mathbf{R}}+\sigma_d^2\mathbf{I}\right)\widehat{\mathbf{w}}_\theta$$

Here's the result comparing the no mismatch, mismatch due to perturbed array and diagonally loaded cases for the MVDR beamformer for the situation in the question. There is an improvement in the signal power estimation, at the expense of increased background noise power.


Method1:
w = w([(end/2+1):end 1:end/2]);

Method#2:
w = circshift(w, ncol/2);


Yes, you could hardcode the values, do an FFT of your signal, and multiply by the hardcoded values, then get the amplitude in the frequency domain from Parseval's theorem (example).

You could also implement a digital filter and do it in the time domain.  Here's Matlab code and my translation to Python.

Note that the bilinear transform causes the digital filter to roll-off more quickly than the analog filter, introducing measurement errors at common sample rates. A simple fix is to upsample the data first, so the discrepancy is minimized. Other filter design approaches here: http://dsp.stackexchange.com/q/36077/29


How to calculate variance of noise on a noisy image (Noise+Signal)

I would assume that I1 and I2 are of the same scene/same lighting conditions/same camera position.

In the ideal world images I1 and I2 would have been related by a constant scale factor (+ random noise) so in this case the way to go would be to determine an optimal scale factor that maps I1 into I2 and then estimating noise variance from the difference. This variance would corresponde to the sum of noise variances of your 2 images, but under the assumption that I1 is significantly less noisy this would give a descent estimate

In the real world there is a thing called camera response curve that maps scene irradiance to pixel intensity values and it's non linear. In fact there can be many more things involved, like gamma correction, custom postprocessing that is done by your particular camera imaging pipeline. Undoing all this mess to get back to linear situation can be an incredible headache. There are software packages and methods for camera response calibration/estimation of other pipeline parameters so in principle it's possible to calibrate the camera prior to taking the images, then undo the camera response curve/other processing making everything linear and after that you can estimate the noise just as in 1). Check camera response curve on google for more details

Applying procedure 1) in the uncalibrated situation gives you a solid overestimate of your noise variance so maybe this would be OK for your purpose.

I had situations where I had to estimate noise level of an image given just a single image. Since typical real world images contain a lot of low frequency regions in practice it's sometimes possible to simply process your image with some smoothing filter (linear or for example median) and then analyse the difference between smoothed image and the original one. For example if at least half of your image correspondes to low frequency regions MAD estimate taken from the difference of your image and smoothed version of your image would give you a usable estimate of noise variance.

Under the assumption of unknown camera response curve you can take a statistical approach to this problem and try to estimate a model with more parameters than just a single scale factor that maps intensities of I1 to I2. Unfortunately I am not ready to give you the exact parametric form of such a model, since I don't see your data.

Hope this helps!


Example exercise - Shannon's law on PCM-digitized signal

I know that for both of those questions I need to use Shannon's law, C = B*log_2(1+SNR) with SNRdb = 10*log10(SNR).

I'm not quite sure that's true. For the first question,

  For SNRdb == 30dB, what is the number of uniform quantum levels needed?

quantization theory is pretty sufficient: In every digitally represented signal, there's so called quantization noise; it's nothing but the power in the difference between analog and quantized signal. For example, every additional bit in the ADC halves the size of the quantization step, reducing the average voltage error by one half, leading to four times less noise power (or, equivalently, 6dB better quantization noise).

So, this question, if you ask me, is: "When you need to achieve quantization noise better than 30dB, how large can your quantization levels be?".

Solution is relatively simple: Let $P_s$ be your signal power. Your Quantization-limited SNR is then

\begin{array}
\mathrm{SNR}_Q &= \frac{P_s}{P_\text{quantization noise}}
\end{array}

Without loss of generality, we set $P_s=1$. Hence, with $q$ the number of quantization steps and 1 the overall ADC range,

\begin{array}\\
{\mathrm{SNR}_Q}^{-1} &= {P_\text{quantization noise}}\\
&= E\left[\left(x_\text{digital} - x_\text{real}\right)^2\right]\\
&= \text{var}({e_Q})\text{, }e_Q\text{ being the quantization noise Random Variable (RV),}\\
&= \int\limits_{-\frac1{2q}}^\frac1{2q}{\left({e_Q}^2q \right)d{e_Q}}\text{, since variance of a uniform RV over }[-\frac1{2q}; \frac1{2q}]\text{ is known,}\\
&= q \left[ \frac13 {e_Q}^3 \right]_{-\frac1{2q}}^\frac1{2q}+C\\
&= \frac q{3}\left(\left(\frac1{2q}\right)^3 -\left(-\frac1{2q}\right)^3  \right)+0\text{, }C\text{ being zero because of plausability}\\
&= \frac q3\left(\frac1{8q^3}+\frac1{8q^3}\right)\\
&= \frac q3\frac1{4q^3}\\
&= \frac1{12} q^{-2}\\
&\overset!\le {30\text{dB}}^{-1}\\
&=\frac1{1000}\\
&\iff\\
q^2&\ge \frac{1000}{12}\\
&\iff\\
q&\ge \sqrt{\frac{500}6}\\
q&\ge 9.1\text{something}\\
q&\ge \left\lceil{9.1}\right\rceil \\
&=10
\end{array}

assuming that the ADC is perfectly linear and thus has absolutely uniform quantization noise for a signal that has uniformly distributed power over the whole ADC range, and thus, the error is uniformly distributed over the quantization step of size $\frac1q$.

That was fun.

  b. What is the minimum data transfer rate needed for the transmission of the digitized-with-PCM signal?

Ha, Shannon, here we come.

PCM means we can't trick around with undersampling. So, our signal's bandwidth is 3kHz, meaning we need to sample with 6kHz. Every sample has information $I=-\log_2(P(x))$, $P(x)$ being the probability that value x is represented in the sample. Still assuming the signal is uniformly distributed (the assignment did not say it was a sum of sinusoids or similar!), $P$ is constant, $P=\frac1q$, hence $H(X)[\text{bit}]= -\log_2{\frac1{10}}= \frac{-\log_{10}{\frac1{10}}}{\log_{10} 2}=\frac1{0.3\dots}\approx 3.32\text b$.

So the bandwidth (bitrate) we need to transport is

\begin{align}
r_b &= 3.32\text b \cdot 6 \text{kHz}\\
&= 19.93\frac{\text{kb}}{\text s}\text{ .}
\end{align}

So with approximately 20kb/s, you can do it.

Those are very good questions. I'd like to ask you to post your own answer if you get an official solution, of course only after getting an O.K. from the professor/lecturer/TA to share his intellectual property (maybe pointing him to this would help). (I of course still like if you accept and upvote this one. It could also certainly use some improvement/corrections!)


explanation of correlation of stationary stochastic processes

You are mixing up two different notions.

Your random process is a collection of random variables $\{X(t)\colon -\infty < t < \infty\}$, one random variable for each time instant. The autocorrelation function of the process is
$$R_X(t, \hat{t}) = E[X(t)X(\hat{t})], -\infty < t, \hat{t} < \infty$$
and is a two-dimensional function (meaning it has two arguments,
and if you try to visualize it, you should imagine a surface above
(and possibly
below) the plane with coordinate axes $t$ and $\hat{t}$. Definitely above
because $R_X(t,t) = E[(X(t))^2]\geq 0$ for all $t$ (and processes for
which $R_X(t,t) = 0$ for all $t$ are kinda boring).

For a stationary process, $R(t,\hat{t})$ depends only on the
difference $\hat{t}-t$ of the two arguments, and it is common
practice to write it as a function of one variable $\tau = \hat{t}-t$.
Thus,
$$\mathcal R(\tau) = R(t,\hat{t}) = E[X(t)X(t+\tau)].$$

The other notion is that if the experiment is performed, then
each random variable takes on a numerical value. Thus, $X(t)$ takes
on a specific numerical value that we name as $x(t)$ while
$X(\hat{t})$ takes on some other numerical value that we name
as $x(\hat{t})$ and so on. Thus, we have a real-valued function
(or signal) $x(t)$ (nothing random about that!) for which the
autocorrelation function is defined as
$$R_x(\tau) = \int_{-\infty}^\infty x(t)x(t+\tau)\,\mathrm dt.$$
Note that $x(t)$ is called a sample path or realization
of the random process: it is what you observe on your oscilloscope
if you lift up your head from the math a little bit
and look at the real world.

It is every DSPer's fondest hope that $R_x(\tau)$ is the same function
as $\mathcal R(\tau)$. Unfortunately, nothing in what has been
said above justifies this hope. But, since people like to assume
that it is indeed true that $R_x(\tau) = \mathcal R_X(\tau)$, mathematicians have come
up with a notion called ergodicity which (for our purposes) says that
for an ergodic process, the ensemble means (expectations)
equal the corresponding time-averages.  As a practical matter,
all model-builders assume that the processes that they want
to model are indeed both stationary and ergodic,
and decide on what the mean and
the autocorrelation function etc should be by making observations of
actual sample paths etc.

For more than what you probably want to know about random processes,
see this answer of mine, and
for an example of a (boring) stationary random process that is
not ergodic, see this other answer of mine.


Fourier transform with periodicity at the harmonic frequency

If you have two periodic functions $f(t)$ and $g(t)$, where $f(t)$ has fundamental frequency $f_0$ and $g(t)$ has fundamental frequency $2f_0$, then their sum $h(t)=f(t)+g(t)$ is periodic with fundamental frequency $f_0$. Without any further knowledge about $f(t)$ or $g(t)$ there is no way to separate the two, because if the only knowledge you have about $f(t)$ is that it is periodic with fundamental frequency $f_0$, then there is no way to distinguish it from the sum function $h(t)$, which has the same periodicity.


Why doesn't JPEG use 1D DCT for image compression?

If JPEG used a 1d DCT it would only be eliminating the redundancy in a single direction. You could compress the image with a 1d DCT then turn it sideways and compress it again. Since the DCT is separable doing this double 1d compression procedure would be (roughly) the same as just using a 2d DCT transform to begin with.


Why can you use the one-sided laplace transform to solve differential equation describing a causal LTI-system?

If it is not equal to 0 for t<0 how can we know we can take the one-sided laplace transform of the LHS?

Unfortunately, the way the problem is framed it is very difficult to do in a rigorous manner.  You can, by remembering that $0^-$ is shorthand for $-\epsilon$ as $\epsilon \to 0$, and taking limits everywhere, and having expressions where you're looking at intervals of $t$ from $-\infty$ to $-2\epsilon$ and other maddening things.
The easy way to do it is to separate the problem: solve the problem for the $x(t)$ that's given.  Then augment the problem by finding an $x(t)$ that is zero everywhere except in $0^- < t < 0^+$ that will result in $y(0^+) = 2$ (note the $0^+$ instead of $0^-$) and $\dot y(0^+) = 1$.  This involves setting $x(t)$ to a linear combination of $\delta(t)$ and $\delta^2(t)$ (as if $\delta(t)$ weren't wacky enough, $\delta^2(t) = d/dt\ \delta(t)$).
Then add the two solutions together.  It may not be what your prof had in mind, but it takes the least amount of hand-waving, and if you really want rigor you can break $\epsilon$ and put it to work.
By doing the above trick of expressing $x(t)$ as a linear combination of the specified $x(t)$ plus whatever sum of the Dirac impulse and it's derivatives, you cast the problem into one where the system is linear -- you're just using physically impossible values for $x(t)$.


Error correcting code? - using RGB space with a lookup table

This is a bit of a "I hope I understand your question right" answer, but seeing I have to leave for the day:
Randomizing the LSBs of all the three color channels would obviously reduce the amount of information a single pixel could carry would reduce the amount of information carried in that color from 24 to 21 bits.
Picking one randomly and randomly varying that has an entropy $1+ \log_2 3$ bits, so you reduce the information contained in that pixel's color from 24 to ca 21.415 bit. This doesn't look like it seriously introduces any privacy there.
It will be impossible to find a 24-bit code that has a rate higher than 21.415/24 that can correct 3 bit errors, so, honestly, you usually simply wouldn't need to use error correction at all.
Still, I think they're doing that for a reason – probably it breaks an existing system of identification. There's no "provable" privacy to be won here, just an "against a single implementation" privacy.

Assuming the anti-aliasing blends linearly

Not the case, usually; you use low-pass filters that are more than a linear interpolation to get acceptably-looking edges.

How to create an ordering of points in 8bit RGB space, such that the minimum number of points are colinear without knowing the final total number of points?

Oh, that's a very interesting problem!
However, we do know from linear algebra how to find linearly independent vectors in a vectors space (e.g. simply choosing canonical base vectors that are all zero but for one entry, or taking a set of random vectors and doing Gram-Schmidt on them).
In this case, our vector space is 3-dimensional. In a three-dimensional space, you get exactly 3 bases in each orthogonal base system, soooo, that's three code points that are orthogonal only! That's very few.
Now, if you add two of these code points, you get another code point, which isn't colinear with any of the existing code points, and since 3-choose-2 is 3, you get three new code points, for a total of 6. You can also add up all three original bases and get one new. OK, so 7 code points.
You can iteratively repeat that, but you have to be careful: You will find multiple ways to calculate the same code point later on, or even multiples of that, and at some point, your elements' coordinates can get larger than your 8-bit integers can hold. But since essentially computers are fast, and the number of things you'll get will be still quite low, you can probably really just brute force your way out of this reasonably by going through all possible pairs of points you have, adding them, and checking that the result is not a multiple of any existing point.


How to detect impulse like peaks - Matlab

Try this:

Convert all values to dB (looks like you've already done this.)
Run the resulting data through a high pass filter with a cut off of 1/4 of your sampling rate.
Select a threshold (looks like 10dB would be good.)  Anything above the threshold is peak.

You'll probably want to use a FIR filter for the highpass - the signal will be delayed depending on the length of your filter, but you can account for that easily since it will be an integer number of samples.

I've been thinking of making an acoustic echogram, and this is the method I thought up to get rid of clutter and show just the sharp reflections.


detect and distinguish driving maneuvers via smartphone sensors

I want to distinguish between harsh acceleration and normal acceleration as well as harsh braking and normal braking (or even coasting).

On the basis of accelerometer readings, you could place thresholds on acceleration regions or look for outliers based on a profile of "normal" driving.

A typical crash deceleration involves tens of gs (fig 7) and airbags deploy after -1.0g for more than 10ms (page 6).

Therefore, anything close to deviations of 1.0g is already looking suspicious in terms of "feelings" to be expected in normal driving conditions. (Of course, you have to make sure that the phone is fixed firmly on the car). As a side note, typical airliner flight does not seem to go too far above 1.0g (of course, it is more during take off and breaking at landing).

The other approach to take would be to carry your device along with you for a large number of "normal" trips and then process the collected waveforms to establish those "acceptable" limits yourself. This could be seen as Anomaly Detection. In other words, you would be deriving a set of features from your waveforms (modal mean of accelleration over a period of time, variance, median, typical waveforms involved with steering or breaking and accellerating, or other) and then evaluate how "close" a given driving session is to those metrics of normality.

  How can I distinguish harsh acceleration and normal acceleration as well as harsh braking and normal braking (or even coasting) in real-time?

Breaking and acceleration are covered above. By coasting, I am assuming that you mean allowing the car to roll in neutral. This is a rather difficult problem to be solved just by looking at a single accellerometer's output. It is very likely that it would have to be augmented by GPS or gyroscope sensors as well (most phones have all three types of sensors today).

To detect coasting, you would basically have to monitor very small changes in the attitude of the chassis. In constant velocity movement, the accelleration is zero and the engine is producing just enough power to counteract the friction of the tyres to the road and the chassis to the air. BUT! Part of the engine's torque is compressing (slightly) the suspension system. Therefore, a car under power tends to "lift its nose" (in side-view, the wheels want to turn in one direction, therefore the chassis tends to turn in the opposite direction. Try to notice the horizon in this video.). Of course, the road is not the track and for every torque setting, it is possible to find just the right downslope that would make the car appear to be coasting but it would still have a gear engaged. Therefore, inclination, either measured by gyro or GPS derived slope would be invaluable in determining coasting.

Hope this helps.


What's the model for a DAC followed by a sampler?

What kind of DAC is that? Standard instrumentation DAC (modeled by a ZOH, with a negligible settling time and no reconstruction filter) or something fancier?

What's happening could be described in the following terms: in the short analog land between the DAC and the rest of the circuit, you see your signal with mirror images of its spectrum above Nyquist frequency (amplitude decreasing with frequency, as your ZOH is equivalent to a multiplication of the spectrum by a sinc function). The switched capacitor circuit resamples this signal; but since its sample rate is the same as that of the original signal, there is no aliasing - the mirror images created by the DAC/ZOH will line up perfectly with the original harmonics.


Relative Computational Complexity for Channelizer Implementations

For sufficiently long filters, it's usually more efficient to implement the individual filters as FFT-based fast convolution, indeed!

Now, whether something is faster than something else really doesn't depend on the FLOPs alone, but also on the kind of operation, the memory / bandwidth requirements and so on.

In fact, this question was raised in the context of GNU Radio's PFB channelizer; luckily, Tom Rondeau wrote a short article on his findings whether direct FIRs or FFT filters in the individual branches are beneficial in the environment of general purpose computers running his PFB implementation.

From that article:

Now, for GPU-based channelizers, this might be a bit different. Not an expert here, but I'd like to point you to Mr. Jan Krämer's pretty nice presentation from FOSDEM17's SDR Devroom on GPU-Enabled Polyphase Filterbanks("Everyday I'm shuffling") (slides here):

The talk highlights that the same number of ops, depending on your hardware, and memory layout, mean a very different throughput.


sound source localization - calculation of sound intensity vector angles

So this technique approximates the velocity of the particles using the pressure gradient at a time instance. This approach is only correct for microphones that are spaced closely together. The max frequency it can measure correctly with the spacer we're using is around 4.8kHz (stated in datasheet), so all angles corresponding to higher frequencies gave a wrong answer. This caused the 180 degree shifted angles.
A great explanation/visualisation from Siemens can be found here


Convolution theorem for cross-correlation

To stick to proper notations, if $\mathcal F \left\{x[n]\right\}=X(\omega)$ and $\mathcal F \left\{h[n]\right\}=H(\omega)$, what you mean is the Fourier transform of a convolution gives the product of the transforms of each term:

$$
\mathcal F \left\{x[n] * h[n]\right\} = \mathcal F\left\{\sum_{k=0}^{\infty}h[k] x[n-k]\right\} =X(\omega)\cdot H(\omega)
$$

As already mentioned, there is a similar formulation for correlation, the correlation theorem. And again, the Fourier transform of the correlation:
$$
\mathcal F\left\{\mathrm{corr}\left(x[n],h[n]\right)\right\}=\mathcal F\left\{r_{xh}[n]\right\}= \mathcal F\left\{\sum_{k=0}^{\infty}h[k] x[n+k] \right\}= X(\omega)\cdot H^*(\omega)
$$

The product $X(\omega) H^*(\omega)$ is called cross-energy density  spectrum and noted $S_{xh}(\omega)$, and this is just $\mathcal F\left\{r_{xh}[n]\right\}$. The $^*$ sign is for conjugate.

Note that your convolution/correlation are on discrete-time signals $x[n]$ and $h[n]$. So $X(\omega)$ and $H(\omega)$ are transforms of $x[n]$ and $h[n]$ respectively, and not of continuous-time signals $x(t)$ and $h(t)$.


Given a causal LTI system with a transfer function, determine if the system is an all-pass filter

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

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


How can I calculate the entropy of a signal that's not independent from itself?

I think that the answer is that there is no correct answer.  Whilst this is a simple question, the answer is kinda open.

You're asking about the entropy of correlated data with respect to random number generation.  The key reference for this is NIST Special Publication 800-90B, Recommendation for the Entropy Sources Used for Random Bit Generation. Look to sections 5 and 6. You be faced with two questions; Is my data correlated (non IID) and if so, what is it's entropy rate? You have already realised that the standard $ -K\sum_{i=1}^n {p_i \log p_i} $ isn't appropriate.

Is your data correlated?

This is the first step, and NIST suggests permutation testing, as below:-

There are 11 tests NIST recommends, all of which are somewhat complex.  To help (not) they have release some code on GitHub that implement these tests.  I'm unsure of how official this is, but see the well received How to interpret the entropy results for a NIST test file as to their soundness. An immediately apparent problem for implementation is that as $ T, T_i \in \mathbb{R} $, the equality test in Fig.4, section 2.2.2 becomes difficult without bounds. In a nutshell, the code doesn't work well enough for field use.

What is the entropy?

If you manage to satisfy yourself that the data is non IID, NIST recommends taking the lowest value of another 10 complex statistical tests as the min. entropy value.  Similarly, the Python and C++ code  has problems outputting a coherent and credible entropy estimate. Perhaps you can implement the tests yourself in another fashion.

What I do.

I too build random number generators and I use similar permutation tests, but focused solely on a compression test. Compression exploits correlation and eliminates redundancy.  Compression is also used in both the NIST non IID confirmation and entropy estimation tests above.

One can calculate what I have termed a correlation factor, $ c_f $, explored in What exactly does compression say about correlation of data?  If $ c_f \simeq 1 $ then the data is IID, and non IID if  $ c_f > 1 $. I then use the non permuted compressed file size and divide twice by two to get a conservative entropy rate.  This method seems to be as good as any, and more robust that most. I estimate about 1.6 bits / byte from a 24V Zener diode at 10kSa/s. That's 16kbps of cryptographic strength entropy.

Notes.

I now use fp8 and paq8px compressors as they are about the best I can find that will compile for me. They'll compress IID data to within 1% of the theoretical Shannon limit.
Divide by two as a security precaution to allow for improvements in compression algorithms, whilst realising that these improvements are asymptotic with diminishing returns.
Divide by two again just because.


Recognition of household objects ( image processing, pattern recognition ) using wide lens?

I share your assumption about depth being useless here.

The approach #1 based on point detectors seems also useless, because there is probably a difference in scale between your reference (learning) images and the representation of objects in the real pictures that is so big that point detection becomes useless. Information is not present at the same level.

So, you're left with approach #2, which sounds reasonable since it's also more-or-less the state-of-the-art in object detection benchmarks.
You may improve the performance of your system by using a training set of pictures that have the size of the objects in the real (test) images.

Also, you will use only the descriptor part of the existing keypoints and apply it in sliding windows. Since oyu have a very poor resolution on objects, I believe that SIFT-like descriptors (e.g., HoG) will have poor performance, while binary descriptors (LBPs, BRISK, FREAK, ORB) will probably behave better (they are less dependant on pixel accurate resolution).

Also, if your objects are very different, you may try a simple image correlation (using again learning images downscaled to the actual size of the objects in the test pictures).


Proof for determining Fourier coefficients

For the periodic signal $x(t)$ with fundamental period $T=1/F_0$, proceeding from where you are:

$$
\Rightarrow \sum^{+\infty}_{k\ =\ -\infty} a_k \left[\int^{T}_{0} e^{j(k-n)\omega_0t}dt\right] = \sum^{+\infty}_{k\ =\ -\infty} a_k \left[\frac{e^{j(k-n)\omega_0t}}{j(k-n)\omega_0}\right]^{T}_{0}
$$

For $k \neq n$ you got $0$ (you can prove that)
And for $k=n$ you see that $\int^{T}_{0}dt = T$

You therefore have:

$$\int^{T}_{0} x(t)e^{-jn\omega_0t} dt= a_kT = a_nT$$

You can proceed from there.

EDIT:

We do not assume the existence of the two possible solutions. The right-hand side of the equation evaluates to that. It can be written in compact form using the unit impulse which has two distinct values depending on the index. (see the Kronecker delta in digital signal processing). Think $\delta_{kn}$.
We still need the summation, but we can start with the integral. It's only that the integral is zero for all values of $k$ except at $k=n$ where it is $T$. So you remain with one term in your summation.

All in all, your initial equation equals: $\displaystyle \sum^{+\infty}_{k\ =\ -\infty} a_kT\delta_{kn}$


Locally Adaptive Thresholding - How to Detect if Image Block is Uniform

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


Evaluating the spectral density of generated noise through the autocovariance

The autocovariance and spectral density are Fourier transform pairs. Therefore, the autocovariance function and spectral density are equivalent.

Have a look to the definitions and examples for white noise in this presentation.


What is the source of the vibrations apparent the FFT

First, I find the 3D waterfall style plots pretty hard to interpret. My preferred method is a 2D color map (e.g. pcolor() in matlab).  It's a lot easier to see your frequencies.

To your first question, it's going to be hard to answer that without knowing a lot more about your system.  Vibration at twice gear mesh frequency is certainly possible.  It's also possible that you have a structural mode.  The first thing to do to investigate this is to do a frequency sweep of your drive speed. Look to see whether this higher frequency tracks with the drive speed.  If it is synchronous with the drive speed, then you know it is definitely related to something rotating.

Your second question is essentially impossible to answer without knowing a ton more information about your system.  e.g. what kind of motor is it?  permanent magnet DC? 3 phase AC induction?  How are you controlling it?  Is there any load on the gearbox or is it spinning freely?  You'd probably need to post all of the electric schematics, pictures of the entire setup, and more data.  I doubt it's the kind of question that would be easily answered on this kind of forum.


What Nyquist says is that $$R_p = 2B\,\text{ Bd},$$ where $R_p$ is the number of pulses per second transmitted, $B$ is the available bandwidth, and the units are Bauds or symbols/sec ($\text{Bd}$). Here it is assumed that the pulses being used are sinc pulses; otherwise, the pulse rate will decrease. In this sense, $2B$ is an upper bound on the pulses per second that can be transmitted. It is correct to say that you can transmit two pulses for every Hz of bandwidth available.

Now, each pulse can carry several bits, by means of its amplitude. If you allow, say, two amplitudes (for example, -1 and 1), then each pulse carries one bit and $R_b=R_p$, where $R_b$ is the bit rate. In general, if you allow $M=2^k$ amplitudes, you'll be able to transmit at a bit rate $R_b=kR_p$.

Shannon capacity does constrain the number of bits that you can transmit per second; furthermore, the usual formulation assumes sinc pulses.


How to prove $\int_{-\infty}^{\infty}x^2(t)dt=\int_{-\infty}^{\infty}\hat{x}^2(t)dt$ where $\hat{x}(t)$ is the Hilbert Transformation of $x(t)$

The Hilbert transform of $x(t)$ denoted by $\hat{x}(t)$ has the following Fourier transform
\begin{equation}
 \mathcal{F}(\hat{x}) = H(w)X(w) = -j \operatorname{sgn}(w)X(w)
\end{equation}
because
\begin{equation}
 \hat{x}(t) = \int\limits_{-\infty}^{\infty} \frac{x(u)}{\pi(t-u)} \ du
\end{equation}
Notice that
\begin{equation}
 \vert \mathcal{F}(\hat{x}) \vert^2 =
 \vert -j \operatorname{sgn}(w)X(w) \vert^2
 =
 \vert X(w) \vert^2
\end{equation}
i.e.
\begin{equation}
\int\limits_{-\infty}^{\infty} \vert \mathcal{F}(\hat{x}) \vert^2
\ dw
=
\int\limits_{-\infty}^{\infty}
\vert X(w) \vert^2
\ dw
\end{equation}
Use Parsevals equation which gives you
\begin{equation}
\int\limits_{-\infty}^{\infty} \vert \hat{x}(t) \vert^2
\ dt
=
\int\limits_{-\infty}^{\infty}
\vert x(t) \vert^2
\ dt
\end{equation}
P.S: For real signals, you could indeed remove the magnitudes (absolute values).


What is the difference between natural response and zero input response?

First it's important to realize that many authors use the terms zero-input response and natural response as synonyms. This convention is used in the corresponding wikipedia article, and for instance also in this book. Even Proakis and Manolakis are not entirely clear about it. In the book you quoted you can find the following sentence on page 97:

[...] the output of the system with zero input is called the zero-input response or natural response.

This suggests that the two terms can be used interchangeably. Further down the page, we find the following sentence:

Thus the zero-input response is a characteristic of
the system itself, and it is also known as the natural or free response of the svstem.

Again, this strongly suggests that the authors believe that both terms are equivalent.
However, on the pages you mentioned they appear to make a difference between the two. And the difference is as follows. The zero-input response is the response which is caused by non-zero initial conditions. It only depends on the system properties and on the values of the initial conditions. The zero-input response is zero if the initial conditions are zero.
The natural response is the part of the total response the shape of which is only determined by the poles of the system, and which doesn't depend on the poles of the (transform of the) input signal. The natural response does depend on the input signal in terms of constants but its form is entirely determined by the system's poles. Unlike the zero-input response, the natural response does not vanish for zero initial conditions.
The total response of the system can be written as the following two sums:

zero-input response + zero-state response
natural response + forced response

The zero-state response is the response for zero initial conditions, and the forced response is the part of the response the form of which is determined by the form of the input signal.
I hope this becomes clear in the following example. Let's investigate the following system:
$$y[n]+ay[n-1]=b^nu[n],\qquad y[-1]=c\tag{1}$$
where $u[n]$ is the unit step sequence. The total response can be computed using $\mathcal{Z}$-transform techniques:
$$y[n]=\left[\frac{1}{a+b}b^{n+1}+\left(c-\frac{1}{a+b}\right)(-a)^{n+1}\right]u[n]\tag{2}$$
The zero-input response is the part of the total response that is determined by the initial condition and that does not depend on $b$:
$$y_{ZI}[n]=c(-a)^{n+1}u[n]\tag{3}$$
Obviously, $y_{ZI}[n]=0$ for $c=y[-1]=0$, i.e., for zero initial condition.
The natural response is the part of the total response the shape of which is determined by the system's pole:
$$y_N[n]=\left(c-\frac{1}{a+b}\right)(-a)^{n+1}u[n]\tag{4}$$
Note that it depends on the initial conditions as well as on the input signal (via the constant $b$).
Also note that it is the shape of the zero-state response that depends on the poles of the system as well as on the poles of the input signal transform. All other responses mentioned here only depend on one of the two sets of poles. The shapes of the zero-input response and of the natural response depend only on the system's poles, whereas the shape of the forced response is determined by the poles of the input signal. The expression for $y[n]$ quoted in your question from Proakis and Manolakis is the zero-state response (because the system is initially at rest), and the first sum is the forced response, and the second sum is the natural response. Since the zero-input response is zero in this case, the sum of natural response and forced response (i.e., the total response) equals the zero-state response
In mathematical terms, the natural response is the homogeneous solution of the difference equation, where the constants are determined such that the sum of the particular solution (the forced response) and the homogeneous solution satisfy the given initial condition. Clearly, the zero-input response is also a solution to the homogeneous equation, but the difference with the natural response is that the zero-input response alone satisfies the initial conditions, because it is combined with the zero-state response, which assumes zero initial conditions. On the other hand, the natural response alone does not satisfy the initial conditions. The initial conditions are satisfied only by combining the natural response with the particular solution of the difference equation (the latter being the forced response).
As mentioned above, we can write the total solution as
$$y[n]=y_{ZI}[n]+y_{ZS}[n]$$
(zero-input response plus zero-state response)
and as
$$y[n]=y_N[n]+y_F[n]$$
(natural response plus forced response). For the given example, we have
$$y_{ZI}[-1]=y[-1]$$
i.e., it is $y_{ZI}[n]$ that takes care of the initial condition. That's also why $y_{ZI}[n]=0$ if the initial condition is zero. $y_{ZI}[n]$ must satisfy the homogeneous equation
$$y_{ZI}[n]+ay_{ZI}[n-1]=0,\qquad y_{ZI}[-1]=y[-1]$$
So if $y[-1]=0$, $y_{ZI}[n]=0$ for all $n$. The natural response also satisfies the homogeneous equation, but not with the initial condition $y_N[-1]=y[-1]$. What is satisfied is $y_N[-1]+y_F[-1]=y[-1]$. This is why the natural response is generally non-zero, even for zero initial conditions. And the natural response is the homogeneous solution which we need to combine with the particular solution (forced response) found in the standard way. We usually have no direct means to find the specific particular solution which, when combined with the special homogeneous solution represented by the zero-input response, will give the complete solution of the difference equation. For this we need another homogeneous solution, and this is the natural response.
Again using the above example will hopefully clarify this. For an exponential forcing signal, the standard (and most straight-forward) way to obtain a particular solution is choosing a scaled version of the forcing function:
$$y_p[n]=Ab^n\tag{A1}$$
(for the sake of simplicity I leave out the unit step $u[n]$, assuming that we consider $n\ge 0$, unless we talk about the initial condition). The constant $A$ is determined by plugging $(A1)$ into the difference equation:
$$Ab^n+aAb^{n-1}=b^n$$
giving $A=\frac{b}{a+b}$. The general form of the homogeneous solution is
$$y_h[n]=B(-a)^n\tag{A2}$$
Of course $y_h[n]=0$ (i.e., $B=0$) is one specific solution, but that's not the one we're looking for. We need to determine the constant $B$ in such a way that the sum of the particular and the homogeneous solution satisfies the initial condition:
$$y[-1]=y_p[-1]+y_h[-1]=\frac{A}{b}-\frac{B}{a}$$
From this equation we get
$$B=\frac{a}{a+b}-ay[-1]$$
which shows that the homogeneous solution we need is non-zero if $y[-1]=0$. $y_p[n]$ and $y_h[n]$ found in this way are identical to the forced response and the natural response, respectively, as shown in $(4)$ and - implicitly - in $(2)$.


How do I reconstruct text from an image using only morphological operations?

If you're willing to add/subtract etc. morphologically transformed images, you can count how many signal pixels are in the vicinity of each pixel, and threshold based on that number.

img = imread('http://i.stack.imgur.com/wicpc.png');

n = false(3);n(4) = 1;
s = false(3);s(6) = 1;
w = false(3);w(2) = 1;
e = false(3);e(8) = 1;

%# note that you could convolve with a cross instead
fourNeighbourCount = imerode(img,n) + imerode(img,s) + imerode(img,w) + imerode(img,3) + img;

%# require at least two neighbours
img = fourNeighbourCount > 1;

If you then convolve with a 3x3 mask that has a hole in the middle, you can get something like this:


Is this how the lowpass butterworth analog filter can be implemented?

Yes but there are a lot of details; like merging R4 with a termination impedance and using pre-distortion to compensate for finite Q in the inductors.  I would suggest reading some wikipedia articles:
https://en.wikipedia.org/wiki/Analogue_filter
https://en.wikipedia.org/wiki/Network_synthesis_filters
https://en.wikipedia.org/wiki/Butterworth_filter
https://en.wikipedia.org/wiki/Electronic_filter_topology#Ladder_topologies
To get a feel for you want.
Then you can decide how far you want to go along this path.  Whether to learn how to do it yourself or find some prepacked program.  But, unless you  feel lucky today, proceed with caution and some knowledge.  There are a lot of in's and out's; pitfalls.
I have always found the field fascinating; but that's partially because it is complicated and fundamentally mathematical :)  Your mileage will vary.
BTW: this is a DSP forum not an analogue design forum.


Find response of discrete time LTI system given input and impulse response

Note that the given definition of $x[n]$ is undefined for $n=0$. But even if you changed it by replacing the unit step $u[n]$ by the shifted step $u[n-1]$ to exclude the value $n=0$, then the problem would still be very hard to solve. So it's more than likely that the correct definition of $x[n]$ has no $n$ in the denominator, and the current form is just a consequence of a typo or of a misinterpretation.


Start with what we know:

The DFT works perfectly for a single cycle waveform input. For a single cycle waveform, all harmonics have an integer number of cycles (by the definition of “harmonic”).
“Spectral leakage” occurs when the input is not a single cycle waveform. As such, at least one frequency component is not represented by an integer number of cycles.

One way to view the DFT is as an algorithm that probes all possible sinusoids of the input waveform, determining frequency and phase of each. Fortunately, we have a limited number of sinusoids to probe for, knowing they must be harmonics.

One way to probe for a specific frequency, if you know the phase to expect, is to multiply by a sinusoid of the same frequency you’re looking for, at the same phase. If you take the average value of that point-by-point multiplication, it gives a value proportional to the amplitude of the component you’re looking for. Against any non-matching harmonic, the result averages to zero. Otherwise, it’s proportional to the amplitude of that frequency in the input. This is true with any number of harmonics present in the input.

Here's an example, probing for the first harmonic. There average of the results of the multiplication, averaged, yields exactly half the peak of the target:

And here is an example of the target not being the same frequency as the probe—the result wave has positive and negative points that average to zero, indicating no match:

The shortcoming for this scenario is that you need to know the phase of the harmonic component in advance. However, an alternative is to probe twice, once with a cosine wave and again with a sine wave. Because of the phase relationship they have, between the two results we can determine both amplitude and phase of the input harmonic under test. Please see my article A gentle introduction to the FFT for the finer details. In the DFT, the real part of the complex “bin” value represents its cosine result, and the imaginary part the sine result.

But, if a frequency component in the input is not an integer number of cycles, it isn’t really a sinusoid and isn’t a harmonic—easy to see if you try to repeat it end to end. And it won’t line up with any of the probe frequencies, and the math can’t have the property of either exactly matching or fully cancelling the individual probes. If the input is very close to a single cycles sinusoid, for instance, but having more than a single cycle, the probe-and-average with the expected first harmonic will give a value a bit less then a perfect match, the subsequent harmonic tests won’t full cancel either. That’s the “leakage”, if you want to view it as the missed energy spilling elsewhere.


Your question is really quite broad.

In a nutshell, if you think of the "successive differences", called spectral flux, understood per frequency range you can start to generalize to different kinds of music where maybe the kick drum is not present but the other instrumentation is sufficient to conclude that there is still on onset. One can also better measure the "strength" (loudness) of an onset, which is useful for algorithms processing onsets to conclude something at a higher rhythmic level (meter, phrasing, other groupings)

There's a lot of literature & code out there to try:

Look into MIREX publications and datasets
Try out code on github (https://github.com/CPJKU/SuperFlux, https://github.com/search?utf8=%E2%9C%93&q=onset+detection&type=Repositories&ref=searchresults, https://github.com/adamstark/Beat-Tracking-Evaluation-Toolbox)
Play with Beatseeker (https://www.ableton.com/en/packs/beatseeker/)


How to build a realtime simple audio spectogram?

The Discrete Fourier Transform of a signal is a series of $\mathbb{C}$omplex numbers. When performed in short time segments and depicted in a spectrogram, what is really depicted is a spectrum.

That is, the strength of a particular frequency component, which is a $\mathbb{R}$eal number.

There are various different expressions for the strength of a component of the spectrum such as Amplitude, Power, Power Spectral Density and so on.

The most straightforward one is the Amplitude spectrum. To calculate the Amplitude of a given spectrum component, you need to get the absolute value of the complex number that represents a frequency component in your spectrum.

Programming languages such as C, Python, Julia, MATLAB / Octave and others, have native complex data types. In that case, the output from an fft() kind of function is an array of type complex and each element of that array contains a function (e.g. abs()) to obtain the absolute value of that complex number. If your programming language does not have a complex data type, then the output of your fft() is probably two arrays of real numbers. One for the real part of the complex number and one for the imaginary part. In that case, to obtain the spectrum, all that you have to do is spc[m] = sqrt(real_part[m]^2 + imaginary_part[m]^2) where ^ denotes exponentiation and m is the $m^{th}$ element of your DFT spectrum.

Once you do that, you will have an array of real numbers. These are going to vary wildly, because of the way the DFT sums work and for this reason you usually apply a logarithmic transformation. This transformation keeps large numbers large and makes small numbers large. You would typically apply this with something like spc[m] = log(1 + spc[m]).

At this point, you can depict spc[m] in grayscale by normalising it to 255 (usually) tones of gray with something like rendered_spc[m] = int(round(spc[m]*255.0)). int denotes an integer data type, you are likely to use something like a BYTE or whatever is more convenient for your bitmaps. The other thing you can do is to map the value of spc[m] within a palette of N colours and "paint" the pixels with that colour.

Finally, it would be probably useful to emphasise here that the spectrogram does not depict one "real time FFT" but many short time fourier transforms. So, if you want to depict 4 seconds worth of a histogram, you would have to accumulate a number of FFTs, get their amplitude spectra, transform them logarithmically, map them to colours and then depict them. Think of each FFT amplitude spectrum as one column of your bitmap.

This is a brief outline on rendering spectrograms practically. Depending on the desired accuracy, you might have to add some more complex features, such as "overlapping" DFT segments, interpolation (for smoother bitmaps with less data) and others.

Hope this helps.