If it's a linear phase FIR filter, then the inputs signals will be shifted (delayed) by an amount of group delay at the output.

For a linear phase FIR filter of length $L = 2K+1$ the group delay will be $N = K$ samples. For even length $L = 2K$ FIR filters it will be at $N = (L-1)/2$ ; half-sample position.

For non-linear phase FIR filters, group delay will be dependent on the specific frequency of the applied input signals.  It should be computed from

$$ \tau = - \frac{d \phi(\omega)}{d\omega} $$

where $\phi(\omega)$ is the phase response of the filter.

The matlab function filter(b,a,x) computes all the samples begining up to input signal length. So the first samples are transients.


That depends, but usually rough frequency correction in an OFDM receiver is done by simply looking at which carriers have significant energy, which don't; you can do that, for example, by comparing the sum of all non-null carriers' energies to that of all null-carrier energies, and shifting that around until the ratio maximizes. Pretty simple idea!


Autocorrelation of discrete signal $\sin{2\pi fn}$

well, there's an old trig identity that you learned pre-calculus:

$$ \sin(\alpha + \beta) = \sin(\alpha)\cos(\beta) + \cos(\alpha)\sin(\beta) $$

or, alternatively

$$ \sin(\alpha) \sin(\beta) = \tfrac12 \big(\cos(\alpha-\beta) - \cos(\alpha+\beta)\big) $$

and you will need to use one of those.

but you have a bigger problem because $x[n] = \sin(2\pi f n)$ is a finite power, infinite energy signal.  the definition for autocorrelation is different, i believe for power signals than for energy signals.

the equation for crosscorrelation for power signals is, i believe,

$$ r_{xy}[m] \ = \ \lim_{N \to \infty} \ \ \tfrac{1}{2N+1} \sum_{n=-N}^{N}x[n] \overline{y[n-m]} $$

for autocorrelation of a purely real signal

$$ r_{xx}[m] \ = \ \lim_{N \to \infty} \ \ \tfrac{1}{2N+1} \sum_{n=-N}^{N}x[n] x[n-m] $$


Determine proper coefficients for the transfer function

Your are basically designing a lowpass filter in order to suppress the noise to some extend. For now I will assume that you are dealing with white noise, so a flat frequency spectrum with a time average of zero. But the value you are after is 26, which in the frequency spectrum is just a spike at zero Hertz. So in order for your filter to recover this value a good as possible it should be stable and have unity gain at zero Hertz.

Fixing $a\in(0,1)$ should ensure that the transfer function has a pole inside the unit disk and therefore be stable. In order to ensure the unity gain at zero Hertz you need to evaluate $H(z)$ at zero Hertz. In order to evaluate a discrete time transfer function at a frequency you can use $z=e^{j\,\omega\,T}$, with $\omega$ the frequency is radians per second and $T$ the sample time. So for zero Hertz/radians per second you get $z=1$. So evaluating $H(z)$ at this $z$ value gives

$$
H(1) = \frac{b}{1 - a}.
$$

This is only equal to one (unity) if $b=1-a$. So if you used $a=0.5$ then indeed it should be that $b=0.5$.


Under what condition is the available bandwidth $f_s$ instead of $\frac{f_s}{2}$?

If your signal is real-valued, then it's spectrum is conjugate symmetric. That means, that negative frequencies (or frequencies from $\frac{f_s}{2}$ up to $f_s$) are mirrored. Thus we can always neglect frequencies above Nyquist range.

Although, if your signal is complex valued, then such symmetry won't exist, and frequencies above $\frac{f_s}{2}$ contain extra information. So you can indeed say that the bandwidth of your signal is $f_s$. This somewhat makes sense, because you might represent now both real and imaginary part of your signal on your frequency axis (twice as much space is necessary).

Here is some small example of sawtooth signal ($f_s=105Hz$), with fundamental frequency of $f_0=10 Hz$ and same signal modulated by the complex exponential with $f_{mod}=5Hz$. In such way we get time domain signal with complex values:
 $$x_{mod}(t)=x(t)e^{2\pi if_{mod}t}$$

As you can see, amplitude (phase also) spectrum for real valued signal is symmetrical around the Nyquist frequency (aliasing is also visible whilst aliased harmonics are not covered by genuine). On the contrary same sawtooth modulated by complex sinusoid (and cosinusoid) is no longer symmetrical. Probably you will notice another "feature" of this signal, while taking a closer look at peak frequencies and relating it to $f_{mod}$...

And traditionally, code:

%% Time domain parameters
fs = 105;
dt = 1/fs;
T = 10*pi;
t = 0:dt:T-dt;
N = length(t);
%% Sawtooth signal generation
f0 = 10;
x = sawtooth(2*pi*f0*t);
X = abs(fft(x))/N;
freqs = (0:N-1).*(fs/N);
subplot(211)
plot(freqs, X)
grid on
hold on
plot([fs/2 fs/2], [0 max(X)], '--r')
xlim([0 max(freqs)])
xlabel('Frequency [Hz]')
ylabel('Amplitude')
title('FT of real-valued signal')
l=legend({'sawtooth','$\frac{f_s}{2}$'});
set(l,'Interpreter','Latex');
%% Complex modulated sawtooth signal
f0_cplx = 15;
y = x.*exp(2*1i*pi*f0_cplx*t);
Y = abs(fft(y))/N;
subplot(212)
plot(freqs, Y)
grid on
hold on
plot([fs/2 fs/2], [0 max(Y)], '--r')
xlim([0 max(freqs)])
xlabel('Frequency [Hz]')
ylabel('Amplitude')
title('FT of complex-valued signal')
l=legend({'complex mod-sawtooth','$\frac{f_s}{2}$'});
set(l,'Interpreter','Latex');


Possible to shift image pixel-wise with FFT?

The point that with "pixel-wise" you mean "differently for each pixel" indicates that a frequency-domain transform is a really bad approach for this. The DFT of a single element is that element (scaled, if at all modified), so you can't win with this.

You can use the FFT and similar algorithms to do elegant shifts, but the fact that it's used in light-field reconstruction stems from the fact that there's a certain, mathematically elegant model of what should happen to many of those pixels that relates nicely to frequency domain.

If you just want to shift/interpolate/distort your image pixel-wise, the normal distortion matrices are the way to go – ask google about that.

For example, OpenCV has a very short tutorial on camera calibration; maybe the theory contained in there is something for you?


How can I implement image widowing in my image in Java?

Consider the large image, and the region of interest (red):

If you will simply take the Fourier transform of the 8x8 region, the transform is assuming that the image is periodic like this:

The Fourier transform will see artificial directional features across the replicates of the 8x8 region. But if you pad the 8x8 image with zeros to 32x32, these artifacts disappear and the frequency resolution of the discrete Fourier transform is enhanced as a welcome side effect as well:

However, now a new problem is revealed: There are artificial directional features along the boundary lines of the 8x8 region, like in the upper half of the left boundary and in the left half of the upper boundary. A solution is to multiply the 8x8 image with a circularly symmetric window function:

The drawback is that the corners of the image are completely ignored. If you think it is useful, you could extend your region of interest so that you would not null the corners of that 8x8 region.

To test the performance of your dominant direction finding algorithm starting with the 8x8 input image, you can use test images. Start with a flat white 8x8 image. This should not give a dominant direction stronger than some level of approximation error. If it does, make sure that your window function has circular symmetry and that you have normalized fairly the calculation of all direction-specific integrals over the squared magnitude of the Fourier transform. Then try 8x8 images with directionality generated by setting the pixels to values $\cos(ax + by + c)$, where $\text{atan2}(a, b)$ defines the directionality, $x$ and $y$ are the pixel coordinates and $c$ should not affect the analysis result.

If your algorithm passes these tests it should be fine!


Power Spectrum Estimation of three sinusoids in white noise

The standard Periodogram uses all the data and computes its power spectrum spectrum estimatation at once; therefore it provides more spectral resolution but less estimation reliability (larger estimation variance).

The Welch's modification to the standard Periodogram is the concept of dividing the signal into shorter blocks and averaging the computed per block Periodograms, and therefore it has less spectral resolution (due to shorter blocks) but provides smaller estimation variance, due to averaging.


Amplitude Spectrum, Nyquist Frequency, mixed/min/max wavelets

Note that moving between minimum phase, maximum phase, and mixed phase involves basically two operation in the time domain:

shifting
time-reversal

If $X(e^{j\omega})$ is the Fourier transform of the sequence $x(n)$ then the transforms of the shifted and the time-reversed sequences are given by

$$x(n-k)\Longleftrightarrow e^{-jk\omega}X(e^{j\omega})\\
x(-n)\Longleftrightarrow X(e^{-j\omega})$$

Since $|e^{-jk\omega}|=1$, the magnitudes of the spectra of $x(n)$ and $x(n-k)$ are identical. Furthermore, if $x(n)$ is real-valued then $X(e^{-j\omega})=X^*(e^{j\omega})$ and consequently the spectra of $x(n)$ and $x(-n)$ also have the same magnitudes.


can $e^{j\phi_0}$ be incorporated into the sine and cosine terms?

Your problem might be the unmentioned fact that the given frequency response $(4)$ in Carlson's book is only valid for positive frequencies. The complete definition should be
$$H(f)=\begin{cases}Ae^{j\phi_0}e^{-j2\pi ft_g},&f>0\\Ae^{-j\phi_0}e^{-j2\pi ft_g},&f<0\end{cases}\tag{1}$$
Now we have $H(f)=H^*(-f)$, which is necessary for the corresponding system (channel) to be real-valued.
If you split the system into two subsystems $H_1(f)=Ae^{j\textrm{sgn}(f)\phi_0}$ and $H_2(f)=e^{-j2\pi ft_g}$, where $H_1(f)$ is a (scaled) phase shifter, and $H_2(f)$ is a pure delay, it's easy to see the impact of the total system on the input signal.
Consider the input signal $x_1(t)\cos(\omega_ct)$. The output of the scaled phase shifter is $Ax_1(t)\cos(\omega_ct+\phi_0)$, if the signal $x_1(t)$ is a lowpass signal with a maximum frequency smaller than the carrier frequency $\omega_c$. Delaying that signal by $t_g$ gives
$$y_1(t)=Ax_1(t-t_g)\cos[\omega_c(t-t_g)+\phi_0]\tag{2}$$
The result for the signal $x_2(t)\sin(\omega_ct)$ follows in a completely analogous way.
In sum, the two necessary conditions for the given formula to be true are

$H(f)$ corresponds to a real-valued system, i.e., $H(f)=H^*(-f)$, and
$x_1(t)$ and $x_2(t)$ are lowpass signals with maximum frequencies smaller than the carrier frequency $\omega_c$.


Its hard to say without knowing exactly how Matlab's randn is processed from the RNG you're using and how Matlab uses the seed.

The most reasonable way is  to save pre-calculated noise data or just deal with a different random number generator (Boost among others has a normal random number generator, and you can pull it into Matlab to use via MEX if necessary).


Apply the Hamming window to a sinusoidal signal

The window isn't "sampled" as it isn't a signal, really.  It's a vector that you construct, and it doesn't matter how you construct it.  The window vector just has to be the same length as your signal vector.


Barcode localization in difficult environments

Other than commercial barcode reading algorithms (many of which fail to read challenging codes), I would like to direct you to this paper which is one of the best academic works in that field in my opinion:

Gallo, Manduchi: Reading Challenging Barcodes with Cameras.

Here is a more recent version:

Gallo, Manduchi: Reading 1D Barcodes with Mobile Phones Using Deformable Templates

Judging by the performance they present it can decode really deformed and low quality barcodes.

Following the localization stage, you will require a robust and fast decoding algorithm. Here is one very strong algorithm from Simsekli and Birdal:

Simsekli, Birdal: A Unified Probabilistic Framework for Robust Decoding of Linear Barcodes.


Calculating RMS/Crest Factor for a Stereo Signal

Conceptually you can think of it as creating one mono signal by concatenating the left and the right channel signals, and then computing the RMS and the Crest factor of that mono signal. With $N$ the number of samples in the current frame, and $x_l[n]$ and $x_r[n]$ the left and right channel signals, respectively, you get the following formulas:

$$\begin{align}x_{\text{rms}}&=\sqrt{\frac{1}{2N}\left[\sum_{n=0}^{N-1}x^2_l[n]+\sum_{n=0}^{N-1}x^2_r[n]\right]}=\sqrt{\frac12\left(x^2_{l,\text{rms}}+x^2_{r,\text{rms}}\right)}\\
C&=\frac{\max\{|x_l[n]|,|x_r[n]|\}}{x_{\text{rms}}}\end{align}$$

where $x_{l,\text{rms}}$ and $x_{r,\text{rms}}$ are the RMS values of the left and right channel, respectively.


Comparing the FFT to numerical integration in Matlab

There are two reasons why you get different results. The first is that you define a vector t between $-10$ and $10$, but the fft command doesn't know that. So the data are interpreted as starting at $0$. This means the first difference is that the FFT computes the transform of a delayed Gaussian impulse, not of an impulse centered at time zero. This is why you get a non-zero imaginary part in the result of the FFT.

The second problem is - as pointed out by nivag in a comment - that the first index of the result of the FFT corresponds to frequency zero, and the negative frequencies appear in the second half of the resulting vector.


A statistical test might be appropriate.  For instance, if you assume the noise floor is Gaussian, then you could compute the standard deviation in any particular octave (or other appropriate sub-band), and look for maxima that are whatever multiple of standard deviations above the local average that you might consider to be statistically significant for your application.


How can I track detected objects from frame to frame

Object detection is relatively a heavy task as you've notice. Detecting the object (in your case human face) in every and each frame would be cumbersome and computationally immense. Therefore, you need to employ an object tracking technique. There are various tracking algorithms, of which, KLT and mean-shift are the two popular ones.
KLT works based on object features (corner,..) but mean-shift (usually) is implemented based on intensity of image.
Personally, I recommend mean-shift for your application, because I think face color would be more unique in a image sequence than general features (corner,edge,..) of the face.
The only thing you might need to add, is verification of the tracker by your object detection algorithm, every like 50 frames to ensure sound functionality of the tracker.

For implementations see:

http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_video/py_meanshift/py_meanshift.html

http://www.pyimagesearch.com/2015/09/14/ball-tracking-with-opencv/


Is there a good method to extract a reference sin/cos from a group of signals?

The best way to accomplish your desired result is to estimate the parameters of the interfering tone(s) and remove it(them) from the signal in the time domain.

To get you started on how to go about this, I recommend you read these two articles of mine:

Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

This needs to be done frame by frame, so there are a whole bunch of implementation details that come into play.  I'm really busy right now, so I will come back later with a fuller explanation.

For the two bin frequency formula, for best results, you will want to define your frame length on an integer + one half number of cycles.  That way you only need to calculate two DFT bins.  The amplitude/phase method will work on just two bins, but you will probably want to do four for that.

Given what you have said, this approach should be very plausible.

I would state that fact as you have 250 samples per cycle.  So suppose you aimed for 2 1/2 cycles, that would mean a frame length of 625 samples.  Any size in that range will put the interfering bin solidly between bin 2 and 3 (zero based indexing) and those are the only DFT bins you need so calculated a full DFT via FFT just for two bins is inefficient.  The longer duration you have, the more bin separation you get between frequencies, but less well the DFT works unless the tones are steady.  All the DFT cares about is sample counts, not what the underlying time scale is.

The proper units for frequency is cycles per frame which corresponds to the bin index or cycles per sample which you have given as 1/250.

Sample code is coming, and so is a further explanation.  You need to think in terms of frames, constructing a copy of your interfering tone in a separate set of buffer because you will need to have overlapped buffers.  Further details also forthcoming.

Still working on my two bin coding.  Hope you are too.

You want to think of your signal as a series of frames.  For each frame you are going to take a reading of the interfering tone.  Some frames will be good reads and some frames will be bad reads, and it is important to have a metric which tells you the quality of each read so you can set a threshold.

You will read the frequency and the amplitude and phase are best done as a cartesian complex number.  You will want to interpolate these values over a span of the good frames to give you parameters to generate a copy of the interfering signal so you can subtract it from your signal.

The reconstruction can be done by overlapping frames with averaging in the time domain as I said before, or a better way is to do a higher degree approximation of the angle in the sinusoidal function.

This method should work quite well for a slowly varying pure tone.  (and possibly the basis of a music compression algorithm?)


Generate 3D image from Multiple 2D Images

You can make a 3D model of the scene (<=> the depth) if the two cameras are looking to the same scene. After that you detect/extract the commun points between the two images and do some triangulation so you can recover the z component.
From your installation as you discribe it you can only enlarge the field of view


Deciding for the correct way to calculate coherence time

This is, as always, a question of what you define something to be. So, context of the different formulas to be key.

For example, the formulas that only differ by a factor of 2 will be taken from two different publications, one understanding bandwidths as double-sided, the other as only single-sided. The definitions of $f_d$ will definitely also differ between these different publications etc.

So, you're not doing your homework 100% honestly: You can't just compare formulas and assume that the same symbols in all of them mean the same thing.

After all, coherence is a term with different meanings for different people: On one hand, a channel is said to be "coherent" over a given time (for a limited bandwidth!) when it doesn't change "too much". Whatever "too much" is for your application. The time within that is the case is called "coherence time".

On the other hand, we often just assign a quantity that we call "coherence time" to a channel, abstracted from whatever you want to do over that channel, based purely on the dynamics of that channel, and those dynamics are captured in the change of the channel over time, which, mathematically, is also a change of phase over time – a freuquency – and hence the term "Doppler Spread" describes that. From that understanding (describing the dynamics of a channel rather than the ability of any specific transceiver to use it) you'll find all the formulas that derive or approximate coherence time from frequency shifts or, rather, bandwidths.

So, the second case applies to your formulas, here. I can't tell you whether that is appropriate for your research – that's up to you as a researcher to figure out (or just: define).

It's (IMHO) very important to realize that this all isn't a deterministic model of what happens to your signal – there's always an underlying, stochastic channel model. For example, you'd model the channel at any instant in time as Rayleigh channel in a tapped delay line model, and then you model the changes to that channel's coefficients with some probability density function. Based on that model, you make a statement about what is "likely" to change (for example, you say "this and that changes with the change following $\mathcal N(0,\sigma^2)$", so let's consider what's within $2\sigma$ as likely)  within a given time (which is the inverse of the bandwidth of the channel). And based on that, you calculate a Doppler spread. So, there's quite a lot of things your channel model defines that contributes to your definition of $T_c$:

A definition of "coherent", typically based on a limit for a (relative) change in some aspect of the channel impulse response
The instantaneous channel model's bandwidth
The stochastic model of how the model changes over time, typically a PDF
An (often arbitrary) definition of what is considered likely to describe any observation over finite time
A way to convert that to a coherence time, typically:

a conversion from some stochastic moment of the PDF to a change-over-time, and an interpretation of that as a frequency
a conversion of that frequency to a time, typically its inverse with some factor depending on normalizations used elsewhere

So, you'll need to pick a channel model first; then things can be derived from, or sensibly defined in accordance to, that model.

There's a lot of good literature out there, but seriously, stochastic, time-variant channels can be a bit daunting at first – take it slowly. Assuming you understand the basics of (deterministic, static) channels: Learn about stochastic, time-invariant channels first, then move on to time-variant channels [1]. The references below are just suggestions. I really can't tell you what works best for you as literature – but if you happen to have a researcher advising you: Researchers have favorite literature that use a notation that they are familiar with, and that is (unlike the wildly collected  formulas you cited) coherent, and they often don't mind discussing things in that literature, so take the chance and ask them for a book on channels that they'd recommend, and read that, and discuss with them.

[1] Proakis: Digital Communications has a chapter "Communication through Fading Multipath Channels" that is the culmination of the channel models it introduces.
[2] Sklar: Digital Communications, as Sklar is traditionally very concerned with a good description of channels; hence, both books introduce channels in the order I recommend above.


what is the amplitude of quantization levels?

Practically speaking an 8-bit linear PCM encoder is implemented inside an 8-bit ADC (Analog to Digital Converter). Such an encoder accepts an analog input within a given refence voltage range and divides this range into 256 equally spaced quantization levels. Some minor deviations may occur based on the quantizer level distributions especially at the first and last levels.

At this stage analog input is converted into a PCM bit stream consisting of 8 bits of quantization codes per each sample.

Now to view these codes ? you can use a binary analyser to see each bit on or off in parallel...

Or if you mean observing the output of a DAC (Digital to Analog Converter) driven by this bit stream then you will be converting these bits back into analog and view it using the oscilloscope. The values will be steppy and be scaled by a reference voltage.


Is this a valid z-transform of $ \frac{1}{n^2}u[n] $?

What you've got is definitely the z-transform of that signal (by definition).  For a more compact expression, notice that the polylogarithm is exactly the infinite sum you've shown (just replace $ z \rightarrow 1/z$ in the Mathworld expression). Unfortunately I don't think there's any nicer way to write it; sometimes there just isn't.


Significance of (pure) sinusoidals in engineering (summary):

1- Periodic signals: Pure sine-waves, represent a useful and simple tool in modeling various physical phenomena, such as standing waves, rotations, alternating currents etc.

2- Fourier analysis: Periodic signals of any engineering sort (not necessarily of pure sine waves) are conveniently decomposed into pure sine waves acording to Fourier series analysis.

3- Linear Systems: strictly speaking complex exponentials $e^{j\omega t}$, are eigenfunctions of LTI systems but sine waves are deeply related. More specifically, the frequency of the sinusoid is preseverd but the amplitude and phase are changed by LTI systems, according to the system's frequency response $H(\omega)$. For instance, for a zero phase system (a symmetric noncausal impulse response) where $H(\omega)$ is real, then the pure sinusoids will also display the eigenfunction property.

4- Nonperiodic Signals: By the extension of Fourier series analysis to integral analysis of continuum frequencies, almost any practical signal is analysed by periodic sine waves.


Transfer function from experimental time domain data

assuming your sets of time-domain data begin with $n=0$, Z transform both the input data, $x[n]$ and the output data, $y[n]$.

$$ X(z) \triangleq \mathcal{Z}\{x[n]\} = \sum\limits_{n=0}^{\infty} x[n] \ z^{-n} $$
$$ Y(z) \triangleq \mathcal{Z}\{y[n]\} = \sum\limits_{n=0}^{\infty} y[n] \ z^{-n} $$

divide the Z transform of the output by the Z transform of the input.  the result of that division is your transfer function.

$$ H(z) \triangleq \frac{Y(z)}{X(z)} = \frac{\sum\limits_{n=0}^{\infty} y[n] \ z^{-n} }{\sum\limits_{n=0}^{\infty} x[n] \ z^{-n} } $$

i cannot guarantee that, if you were to factor numerator and denominator, there will not be identical poles with zeros that could cancel and make the transfer function simpler, but it's a transfer function.


Noise amplification with inverse filtering

The noise amplification would be the ratio of the norms:

$\frac{\text{Norm of noise after inverse filter}}{\text{Norm of noise before inverse filter}} \leq \frac{||\mathbf{A}^{-1}||*||\mathbf{w}||}{||\mathbf{w}||}=||\mathbf{A}^{-1}||$

It seems your reasoning is correct that the norm of the inverse filter is an upper bound on the amplification.


What exactly happens with Spread Spectrum Modulation

Direct-sequence spread spectrum (DSSS) is a technique that is used to generate a modulated signal that occupies more bandwidth than would be implied by its information content alone. A DSSS signal is generated via multiplication of a (typically digitally-modulated) baseband signal by another spreading code waveform. The spreading code waveform is constructed such that it is modulated at a higher rate than the baseband signal; this yields the "spreading" of the resulting signal's spectrum.

To see this, recall the multiplication property of the Fourier transform. Let the baseband signal $x(t)$ be multiplied by a spreading waveform $s(t)$:

$$
x(t)s(t)\Longleftrightarrow X(f) * S(f)
$$

where the $*$ denotes convolution. If $S(f)$ is broadband when compared to $X(f)$, then the result of convolving $X(f)$ with $S(f)$ is that $X(f)$'s spectral content is "smeared" or "spread" across a wider region of the spectrum. The amount of spreading is dependent upon $S(f)$ and therefore the choice of the spreading sequence $s(t)$. So, for example, you could take a 1 kbps BPSK signal, which would normally have a null-to-null bandwidth of 2 kHz, and spread it to instead have a null-to-null bandwidth of your choice: say, 1 MHz.

To answer your most pointed question, you're not really sending a single signal redundantly on multiple carrier frequencies, nor are you dividing the input signal into pieces that get transmitted at multiple carriers. In lieu of a discrete multi-carrier approach, DSSS instead is more of a continuous smearing of the input signal over a large region of spectrum.

If the receiver is also cognizant of the spreading sequence $s(t)$, then it is capable of despreading its observation of the transmission and therefore recovering the original baseband signal. As long as the spreading sequence $s(t)$ does equal zero for any $t$, then in principle, the spreading process is perfectly reversible. In practice, this isn't exactly the case, as there will be some distortion between what is transmitted and what the receiver sees, but it leads to an important point that is often misunderstood about DSSS:

  No benefit in bit error performance, measured in $\frac{E_b}{N_0}$ (the typical metric used for digital modulation performance evaluation) is offered by using direct-sequence spread spectrum. Since the spreading process is reversible, then the principle of reversibility states that there can be no improvement in performance of the optimal receiver for the spread signal when compared to the performance of the optimal receiver for the original, unspread signal.

So if there's no direct bit-error performance benefit, why would you want to do this? There are a few common reasons:

Multiple access: Code division multiple access schemes work using DSSS techniques. CDMA allows multiple transmitted signals to occupy the same bandwidth over the same time period. Transmitters are associated with separate carefully-selected spreading codes; these are also known by the receiver, which uses them to separate the previously-intermingled-in-time-and-frequency signals. CDMA can be made to seem like the greatest thing since sliced bread, but there are some practical issues (most notably the near-far problem) that limit its usefulness. It is, however, used in some cellular networks around the world, for one example.
Interference mitigation: Signals with high bandwidth are resistant to narrowband interferers on the AWGN channel. A hand-waving time-domain explanation: a wideband signal in the frequency domain will have a narrow autocorrelation function in the time domain, and vice versa. The optimal receiver will use a correlator (often implemented with a matched filter) matched to the transmitted spread-spectrum waveform, so its output (which ideally will contain copies of the signal's autocorrelation function at symbol-spaced intervals) will contain narrow, sharp peaks that will be easy to identify against the oppositely-characterized background interference.

A hand-waving frequency-domain explanation: the spreading process takes the original baseband signal's power and spreads it across a wide region of the spectrum. In the communication channel, narrowband interference finds its way into the signal that eventually is observed. During the despreading process, the DSSS signal is compressed back to the original baseband signal's bandwidth. However, at the same time, the despreading process actually spreads the narrowband interference across the same bandwidth as the original DSSS signal.

So effectively, the narrowband interferer has been transformed into a wideband interferer across the entire DSSS signal bandwidth, which as you recall, is much larger than the baseband (information-carrying) signal's bandwidth. Therefore, much of the interference power is pushed out of band with respect to the signal of interest, which yields a higher signal-to-interference power ratio (and better performance).
"Discreetness": No, not discrete, discreet. DSSS is often spoken of as being able to "pull signal out of the noise floor." This refers to the ability of a spread-spectrum receiver to transform a wideband DSSS signal back to the original, much narrower waveform. For a given transmitter power, as the transmitted spread signal bandwidth increases, the average power spectral density over the DSSS bandwidth will decrease, perhaps to the point where it is difficult to observe using spectral analysis at the receiver. This property can be useful for a number of reasons.

One aim might be to decrease the probability of the emitter being detected (and possibly intercepted) by an adversary. Another might be to mitigate interference with other, narrowerband systems that operate in the same region of spectrum. An example of this approach is in ultra wideband communications systems, which have increased in popularity in recent years.


The concept of multiresolution has roots on the observation that in data (eg images), objects and features can be observed and processed at different resolutions. One of the crudest version is to build a pyramid with the original image, and other sub-sampled versions, keeping 1 out of 4, 1 out of 16... pixels.
This yields a crude multiresolution pyramid, which is not a wavelet per se.

Wavelets, in a strict sense, naturally yield a certain type of multiresolution, but satisfy other constraints (esp. invertibility), so they are only a quite useful way of doing multiresolution.

There are many fast ways to compute CWT, that differ from DWT, see for instance

Should I ever pick the continuous wavelet transform over the discrete one? DWT vs CWT vs STFT
Implementing Continuous Wavelet Transform


Best to base it on the highest frequency response you anticipate for the system and implement twice that frequency with the samples you input and the data you read.


How to go about implementing a program to measure the "intensity"/"saturation" of a color in an image?

You can find mapping from RGB space to HSV here:
https://math.stackexchange.com/questions/556341/rgb-to-hsv-color-conversion-algorithm

For other spaces, you can refer to vision modelling textbooks, i.e: Digital video quality, Stefan Winkler


use of complex number in digital signal processing

Any complex array could be represented as a real array as follows:
$$x_R = \begin{bmatrix}
 \Re(x) \\
\Im(x)
\end{bmatrix}$$
and matrices are represented as
$$A_R = \begin{bmatrix}
 \Re(A) & \Im(A) \\
-\Im(A) & \Re(A)
\end{bmatrix}$$
where $\Re(x),\Im(x)$ are vectors containing the real and imaginary part of complex vector $x$, respectively. Note that addition/subtraction is preserved by simply addition
$$x_R + y_R = z_R$$
You can also verify that matrix multiplications in the real domain also preserves the real/complex information in their respective blocks, namely
$$A_RB_R = C_R$$
where $C_R = \begin{bmatrix}
 \Re(C) & \Im(C) \\
-\Im(C) & \Re(C)
\end{bmatrix}$


For small $\Theta$ an approximation would be to simply halve the imaginary value.

You could also use the ratio of the real and imaginary values as the index to a lookup table that would return the half angle value.  You could not store every possible value, of course, so you would have to choose some set of values and then lookup the value that is closest to the actual ratio.


Determine stability of feedback system from open loop transfer function and Nyquist stability criterion gives different results

However when I look at the closed loop transfer function, I would say that this system is unstable for =−1. In this case the transfer function becomes infinity so a bounded input will result in a unbounded (=infinity) output.

This depends on your definition of stability.  $GH = -1$ is called marginally stable because depending on how you look at it, it could be stable, or it could be unstable.
In the Lyapunov sense that @Petrus1904 mentions, it is stable.  But if you take bounded-input bounded-output (BIBO) stability to mean that the input can go on infinitely long but within bounds, then the output can, indeed, go to infinity.  So in the BIBO sense it's unstable.
Hence, "marginally stable".

In my train of thought the point =−2 would again be stable since  will be finite again, however conform the Nyquist stability criterion this point will still be unstable?

An unstable system can have a transfer function that never goes infinite for $s \in j \omega$.  All it requires is that there are no poles on the stability boundary.
I apologize for not taking the time to do so, but I know that I could dummy up a system that has $GH = -2$ and was still stable, because I've designed such systems.  All you need to do is wrap a double integrator with a PID controller: let $G = 1/s^2$ and $H = k_i / s + k_p + \frac{s}{\tau_d s + 1} k_d$ and tune for stability.  The resulting system will have a spot in its open-loop Bode plot which has a phase shift of 180$^\circ$ and a gain greater than one -- to get that gain equal to two you'll just need to jigger the numbers around.
Addendum:  Use the PID transfer function I give above with $k_i = 0.002$, $k_p = 0.02$, $k_d = 0.2$ and $\tau_d = 0.01$ and you should get a system that's stable, with $GH \simeq -2$ at somewhere around $\omega = 0.1 \mathrm{\frac{rad}{sec}}$.
Assuming I got my math right...


You put zeros in between to show delay between the channel taps. I don't know what $M$ is in your code but suppose $M=10$ and say your sample rate is 10 MHz. Then you can interpret mc = [1 zeros(1, M) 0.28 zeros(1, 2.3*M) 0.11];  as you get the first multipath component with zero delay (gain = 1), then 10 zeros later (which is equal to 1 microsecond) you get another multipath component (gain = 0.28), and finally after another 2.8 microseconds you receive the last multipath component (gain = 0.11).
It looks like the "mild" vs "harsh" is just language that the author uses to let you know that mild happens to rarely cause errors where the harsh causes many errors.


it doesn't tell you new information, but because the abs value function is non-analytic (i.e. not all of its derivatives are continuous), and the magnitude-squared is, then the latter can be manipulated mathematically in ways that the former cannot.  one important property of the power spectrum is that it is the Fourier transform of the autocorrelation function in the time domain.

the PSD gives you the energy of each frequency component for a deterministic signal and the expectation value of the energy for each frequency component for a random (or "stochastic" if you like jargon) signal.


How to get the channel impulse response by using PRBS signal with BPSK modulation?

Do not think of your transmitter as BPSK since you are not intending to do data transmission (as you will not be doing data demodulation in the receiver) but simply an upconverter of a bipolar waveform (yes as far as the trasmitter operation goes, this is indeed identical to BPSK but I think that is giving you much confusion). So the transmitter is frequency shifting the baseband signal to an arbitrary carrier frequency.

The receiver simply down-converts the signal back to baseband.

Because of the phase and frequency offsets (in addition to the multipath distortion) the signal as received will no longer exist at baseband as a real signal but will have real and imaginary components. We need to maintain all of these components in order to solve for the channel. Therefore you must do a quadrature down-conversion that would maintain the real and imaginary components (I and Q) of the complex waveform. A very simple example is a channel with an impulse response represented by an impulse at t = 1 sec with a phase of 45°.

$$c(t) = \delta(t-e^{j\pi/4})$$

The received signal is the convolution of our transmitted signal with the channel, so in this case the BPSK transmitted signal would be received as a complex having a 45° phase rotation:

This is a very simple channel; in practice you will have multiple delayed copies at different delays and amplitudes combining in a possibly unrecognizable pattern in the receiver. If the transmitted waveform and the resulting received waveform is known, then you can use inverse convolution techniques such as the Wiener-Hopf equations to solve for the unknown channel, I describe this in detail in the following post:

Compensating Loudspeaker frequency response in an audio signal

(For the case of solving for the channel instead of the channel compensator, swap Tx and Rx as described in that post)

A PBRS (Pseudo-Random Binary Sequence) is an excellent choice for a training sequence as it has spectral content across the entire band. The equalizer can only determine its compensation at frequency locations where energy is present.


Fundamental misunderstanding of windowing

Your description of windowing a 2D image filter kernel, fits into a discrete-time FIR filter design method. Hence I will assume you are mixing windowing as a filter design method vs windowing applied to signals for Short-Time Fourier analysis or windowing applied for block-based processing of long signals. So please make it clear which action you are performing and the context of your intentions.

Then let me briefly describe windowing as a filter design tool.

An ideal discrete-time frequency selective LTI filter will have an infinite length and noncausal impulse response $h_i[n]$. So it shall be truncated and delayed for any practical FIR convolution applications.

Given an ideal filter impulse response $h_i[n]$, you obtain the coefficients of the practical impulse response $h[n]$ of length N samples via

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

Which is the point wise multiplication of the filter $h_i[n]$ with a window sequence $w[n]$ of length N. The main purposes of the window is to first truncate the infinite length ideal impulse response, and then to try to cure the side effects on the resulting filter spectrum caused by the crude truncation.

There are many different types of windows such as box-car,triangular, Hamming, Hanning, Blackman, Kaiser. The box-car is the crude rectangular window without any tapering. Other windows perform varying degrees of tapering at the ends of the impulse response so as to achieve a tradeoff between filter transition bandwith and passpand / stopband peak ripples.

Note that the windows's frequency response $W(e^{j\omega})$ is convolved with the ideal filter frequency response $H_i(e^{j\omega})$ which yields the frequency response of the practical filter;i.e,

$$H(e^{j\omega}) = \frac{1}{2\pi} H_i(e^{j\omega}) \star W(e^{j\omega}) $$.


Why doesn't sampling a periodic continuous-time signal yield a periodic discrete-time signal?

If the ratio between your sampling frequency and the frequency of your signal is irrational, you will not have a periodic discrete signal.

Assuming you have a 1-kHz sine wave and you sample at 3000*sqrt(2) Hz. You will have approximately 4.2 samples per period. However you will not be able to sample the sine wave exactly at the same place. Hence your digital signal will not be periodic.

However, if you sampled the same 1-kHz signal at 4 kHz, you would get a periodic discrete signal. The period would be 4 samples.


Time and frequency domain spreading in CDMA and MC-CDMA?

The spectrum of DS-CDMA without pulse shaping is equivalent to the sinc function. The Matlab code below illustrates the difference between DS and MC-CDMA. For MC you need the IFFT not the FFT function.

clear,clc

N=2;
% BPSK data
d=sign(randn(1,N));

% 7-chip m-sequence code
c=[1 1 1 -1 1 -1 -1];

% DS-CDMA
s1=kron(d,c);

% MC-CDMA
L=sqrt(length(c));
s2=L*ifft(d(1)*c);
s3=L*ifft(d(2)*c);
% etc, better to use a for loop for more symbols


Computing impulse response for system with feedback terms

Your values for $h(0)$, $h(1)$ and $h(2)$ are correct. However, the general formula isn't.
It should be

$$h(n) = \begin{cases}
0, & n<0 \\
\frac{1}{2}, & n=0 \\
\frac{3}{2}\left(\frac{1}{2}\right)^n, & n > 0
\end{cases}$$

A simple way to solve such problems is to use the $\mathcal{Z}$-transform:

$$Y(z)=\frac{1}{2}X(z)(1+z^{-1})+\frac{1}{2}Y(z)z^{-1}$$

which gives

$$H(z)=\frac{Y(z)}{X(z)}=\frac{1}{2}\frac{1+z^{-1}}{1-\frac{1}{2}z^{-1}}=
\frac{1}{2}\frac{1}{1-\frac{1}{2}z^{-1}}+\frac{1}{2}\frac{z^{-1}}{1-\frac{1}{2}z^{-1}}$$

With the transform pair

$$\frac{1}{1-az^{-1}}\Longleftrightarrow a^nu(n),\quad |a|<1$$

and with the property

$$x(n-1)\Longleftrightarrow z^{-1}X(z)$$

you arrive at

$$h(n)=\frac{1}{2}\left(\frac{1}{2}\right)^nu(n)+\frac{1}{2}\left(\frac{1}{2}\right)^{n-1}u(n-1)$$

which corresponds to above result.


Equalizing audio under linux / Raspberry Pi

You can simply implement this as an IIR filter directly in Python using https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfilt.html
Each Room EQ wizard section can be turned into a second order filter section using RBJs Audio EQ Cookbook (See for example https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html or https://www.musicdsp.org/en/latest/Filters/197-rbj-audio-eq-cookbook.html). You want the "peaking EQ". Then just cascade all sections into a single SOS.
Note: If your filter parameters look like https://github.com/mikebrady/shairport-sync/wiki/Digital-Signal-Processing-with-Shairport-Sync I advise some caution: You have a lot of VERY high Q features are 6 Hz/12 Hz (very low freqencies). These could be numerically unstable and you really don't need them. I would eliminate all sections below 30 Hz and just make sure you have decent highpass filter somewhere in the signal chain.


No of samples 4mega samples per second,using 8 bit ADC,V p_p is6V. calculate frequency and peak power

There are sam samples per second and 8 times less pure sine cycles per second. This means that there are sam/8 pure sine cycles per second.

What shannon has to do with that? Probably you mean that Niquitz allows you to reduce sampling, make it 4 times more rarely, as low as sam/4 samples per second. That would sample twice per sine period.


Why is the DFT used for spectrograms rather than the DCT?

Let's define the N-point IDFT $y[n]$ of a signal $Y[f]$ as

$$\begin{align*}
y[n] &= \sum\limits_{f=0}^{N-1} Y[f] e^{j2\pi n\,\frac{f}{N}},& n\in \{0,\dots,N-1 \}\tag{1}
\end{align*}$$
The DCT-II (which is most probably what we're looking at; the others are mathematically useful, but less nice to implement) $\mathbf{Y}[f]$:

$$\begin{align*}
\mathbf{Y}[f] &= \sum\limits_{n=0}^{N-1} y[n] \cos\left(\pi (k+\frac12) \frac nN \right) ,& f\in \{0,\dots,N-1\}\\
&\text{which, Euler says, is}\\
&= \frac12 \sum\limits_{n=0}^{N-1} y[n]\left(e^{j\pi (k+\frac12) \frac nN } + e^{-j\pi (k+\frac12) \frac nN }\right),& f\in \{0,\dots,N-1\}\text.\tag{2}
\end{align*}$$

Now, calculating the inverse discrete Fourier of the discrete cosine transform becomes but a matter of diligence; calculating by inserting $(2)$ into $(1)$, you'll notice that the complex sinusoids are orthogonal w.r.t. the scalar product, ie. the sum over the element-wise products becomes zero for everything but the same exponents.


How does Overlap-Add work for IIR filter?

According to wiki, it seems that overlap add technique usually works for FIR filters. I hasn't read the link you posted and the matrix form for IIR before. But when it comes to buffer for IIR, I have an alternative (hope simpler) method which might works in your case.

  If I want to process a buffer of audio with the IIR filter coming from the LPC analysis, the IR should be infinite, so the only way I can think of using a block convolution is by convolving

The concept is filtering in sections and utilizing the memory in filters. Let's first look at the testing code and usage with following steps,

Dividing the x in to x1 and x2,
Filtering x1 and x2 seprately (with memories in filterszf)
Combining the ouput together.

The results confirm that filtering in section can be the same as x directly.

%% Filter in section
x = randn(10000,1); % psuido signal

x1 = x(1:5000);
x2 = x(5001:end);

b = [2,3];
a = [1,0.2]; % IIR filter coeffs

[y1,zf] = myFilter(b,a,x1); % The zf is the memory in IIR filter!
y2 = myFilter(b,a,x2,zf); % assign the zf when filter the next section.

y = myFilter(b,a,x); % filtering x directly.

isequal(y,[y1;y2]) % Verification

Here is the implementation of myFilter Function. In fact, myFilter imitates(implements) the built-in function fitler in matlab. The implementation architecture is Direct form II, which has advantages over least usage of delay operator and disadvantages over not free for numerical overflow for fixed points.

function [Y, z] = myFilter(b, a, X, z)
    n = max(length(a),length(b));
    a = [a zeros(1,n-length(a))];
    b = [b zeros(1,n-length(b))];
    z(n) = 0;
    b = b / a(1);
    a = a / a(1);
    Y    = zeros(size(X));
    for m = 1:length(Y)
       Y(m) = b(1) * X(m) + z(1);
       for i = 2:n
          z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
       end
    end
    z = z(1:n - 1);
end


One good way is to utilize Bezier curves (or Catmull Rom splines). This way, you could define the control points in the way you like. Positioning of the control points would then allow you to adjust the histogram to the desired output.

I would suggest Catmull Rom splines if you want them to pass through the control points. It gives you better control over the curve behavior. Use 5 such points: 2 for the end points, 1 for the middle, and 1 on each side of the center to give you certain curvature.

If you don't want to use splines, the black curve looks like a Sigmoid function, while the yellow is a third order polynomial.

You can checkout this code to see different implementations.

If you are keen on using Logistic family, then I would suggest the generalized logistic growth curve with horizontal asymptotes at $−1$ and $1$:

$y(t)=-1+\frac{2}{[1+ke^{-\beta t}]^{1/v}}$

where $k > 0$, $β > 0$ and $ν > 0$. The parameter $k$ can be adjusted by a time shift, $β$
corresponds to a scaling of time; $ν$ is the key parameter that determines the shape of the
growth, the asymmetry depending on whether it is less than or greater than $1$. This formula reduces to the standard logistic growth for $ν = 1$, $k = 1$, $β = 2$. The turning point you refer to is commonly known as the inflection point, where the derivative changes sign.

This is explained in detail Ayse Humeyra's article and here. Here is a MATLAB code:

function [y]= general_logistic(t, k, beta, v)
y = -1 + 2./((1+k*exp(-beta.*t)).^(1./v));
end

Finally, in Dino Dini's blog, he explains the normalized tunable sigmoid functions. There, the solution to your problem is addressed. Moreover, a similar question is already posted and answered here. Last but not least, polynomial shaping functions, which would serve the same necessities are discussed in Golan Levin's page.

I hope all these would help.


The phase delay to time delay conversion depends on the frequency of the signal for which that phase delay is associated.

Note that unless the phase delay is given as correctly "un-wrapped", the absolute time delay could be ambiguous by integer multiples of the period of the frequency at which the phase delay was specified or measured.


HINT:

From the definition of the DTFT

$$H(\omega)=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}\tag{1}$$

derive the following facts:

$H(\omega)=H(\omega+2\pi)$
$H(\omega)=H^*(-\omega)$ for real-valued $h[n]$

where $*$ means complex conjugation. Combine these two results to prove the given equation.


Good test for periodicity between signals

Something you could do is calculate the FFT of both signals and define some criterion for periodicity. For example, you pick the highest component of each FFT and compare the two (one from each signal) and thus determine which is more "periodic". You could consider an average of all the FFT components and compare. This would be a rather robust solution, since the FFT calculates the frequency decomposition of a signal.


Does Gamma Correction Interfere with Sub Pixel Motion Estimation?

When one does an operation Wrong / Right aren't strictly defined.
In most cases the questions are:

What's the model?
Is the model reasonable?

When you do sub pixel motion estimation on frames which are highly correlated you can assumes the change in values, even after non linear operation like Gamma Function, is small.

So basically the model assumes that even after Gamma Correction the spatial correlation between 2 frames is kept.

You can do a small test by generating image supposedly in Linear Color Space.
Generate another image as a result of a small movement of the previous image.
Have a look on the spatial correlation between the 2.
Now transform each using Gamma Function and measure the correlation as well.

If their model is reasonable, the correlation in both cases will be similar.


How to assess perceptual volume from FFT

If you equate perceptual volume to loudness, there are international standards.

if you look at
https://www.mathworks.com/matlabcentral/fileexchange/46819-a-weighting-filter-with-matlab-implementation

This is an A-weighted Matlab FFT based implementation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  A-weighting Filter                  %
%              with MATLAB Implementation              %
%                                                      %
% Author: M.Sc. Eng. Hristo Zhivomirov        06/01/14 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xA = filterA(x, fs)

% function: xA = filterA(x, fs)
% x - signal in the time domain
% fs - sampling frequency, Hz
% xA - filtered signal in the time domain

% determine the signal size
sz = size(x);

% represent x as column-vector
x = x(:);

% signal length
xlen = length(x);

% number of unique points
NumUniquePts = ceil((xlen+1)/2);

% FFT
X = fft(x);

% fft is symmetric, throw away second half
X = X(1:NumUniquePts);

% frequency vector with NumUniquePts points
f = (0:NumUniquePts-1)*fs/xlen;

% A-weighting filter coefficients
c1 = 3.5041384e16;
c2 = 20.598997^2;
c3 = 107.65265^2;
c4 = 737.86223^2;
c5 = 12194.217^2;

% evaluate the A-weighting filter in the frequency domain
f = f.^2;
num = c1*f.^4;
den = ((c2+f).^2) .* (c3+f) .* (c4+f) .* ((c5+f).^2);
A = num./den;
A = A(:);

% filtering in the frequency domain
XA = X.*A;

% reconstruct the whole spectrum
if rem(xlen, 2)                     % odd xlen excludes the Nyquist point
    XA = [XA; conj(XA(end:-1:2))];
else                                % even xlen includes the Nyquist point
    XA = [XA; conj(XA(end-1:-1:2))];
end

% IFFT
xA = real(ifft(XA));

% represent the filtered signal in the form of the original one
xA = reshape(xA, sz);

end

The Mathwork will also sell you their Audio toolbox that has many types of weightings.

I find these standards a bit confusing because you don't know if they correct for the effective noise bandwidth of proportional bandwidth filters and if an FFT implementation should take that into account.

I don't vouch for the Matlab code above.  It seems reasonable.

The Wikipedia article:

https://en.wikipedia.org/wiki/A-weighting

also gives the analog transfer functions, which can be converted to digital filters, but actual digital filters seem harder to find.

If your ultimate purpose is to use a standard where a professional signature is required, or one is exposed to liability, I would be inclined to use the Mathwork's product.


Degraded image: Wiener Filter vs Wiener deconvolution

Wiener deconvolution is an approach to solve the deconvolution problem that relies on the filter proposed by Wiener. The equation is the same in denoising and deblurring, except that the filter $G$ (to stick with Wikipedia's notations) that you should use is different.

To make things clear:

denoising consists in the case where the degradation kernel $H$ is the identity,
in deconvolution this kernel is equal to the blur kernel that destroyed your image (motion
blur, lens PSF...).

Hence, in denoising, the $H$ part that appears in the expression of the Wiener filter $G$ disappears (it is equal to 1), while in deconvolution it is equal to the Fourier transform of your blur kernel.

More generally, when using an inverse problem approach for these tasks, i.e. solving $\| A x - y \|^2 + \lambda g(x)$ (with $x$ your solution and $y$ the noisy observation, $g(x)$ some prior on $x$), you can always switch between denoising and deconvolution by choosing $A = Id$.


Machine learning with only positive labels

I would use a novelty detection approach: Use SVMs (one-class) to find a hyperplane around the existing positive samples. Alternatively, you could use GMMs to fit multiple hyper-ellipsoids to enclose the positive examples. Then given a test image, for the case of SVMs, you check whether this falls within the hyperplane or not. For GMMs, you check if it is enclosed in the hyper-ellipsoids. They are both proven to work well in practice.

If you also have some unlabled data in your training set, I would certainly adapt a variant of transfer learning. Maybe you would be able to automatically label the unlabeled data, based on the already learnt samples.


Relationship between "received" signal power and "transmitted" energy per bit?

You are correct that it could be argued that every physical channel will introduce attenuation. In digital communications, however, it is common to define "channel models" that simplify the equations and make the problem easier to understand and to solve.

The most popular model is the "additve white Gaussian noise" (AWGN) model, where the only effect of the channel is to introduce noise. In other words, there is no attenuation. In that model, Sklar's relationship between power and bit energy is correct. This model is very popular because, in many cases, you can ignore the attenuation by assuming an appropriate gain mechanism in the receiver.

Of course, in many other cases the AWGN model is not appropriate. In the wireless channel, for instance, the channel gain is so fundamental that it gets incorporated into the model. If the channel is narrow-band, then the relationship between power and bit energy is similar to what you propose, with $\alpha$ random and possible time-varying.


Compression algorithms specific to complex signals

Complex signals are a special case of multidimensonal signals (where the dimension is two).
A lossy approach tackling compression of multidimensional signals is vector quantization. A very good resource is the book:

"Vector Quantization and Signal Compression",

co-authored by Robert M. Gray.

Vector Qquantization is a classic lossy source coding technique that is used to quantize $M$-dimensional source vectors into a set of $N$-dimensional code vectors which is referred to as the codebook, and $N<M$. A scalar quantizer is a special case with $M = 1$.

A complex signal in general can be seen as two correlated sources. Of course, if you don't want to assume any correlation between the real and imaginary part, then the answer is trivial since the two sources can be compressed independently and separately (using the standard compression techniques of real signals).

If you are interested in information-theoretic aspects of lossy/lossless compression of correlated sources, Wyner-Ziv and Slepian-Wolf coding are the keywords. These schemes are refereed to as distributed source coding schemes. However, since a complex signal can be modeled by two correlated sources, they can be considered here as well.

Finally, for analytic signals, since the imaginary/real part can be perfectly reconstructed from the other one by using the hilbert transform, the compression problem boils down to compression of a single real signal, for which there are plenty of standard schemes.

In response to the updated question:

Loyd algorithm generalized for vector quentization is referred to as LBG Algorithm (abbreviation of the name of the authors of the following paper)

Y. Linde, A. Buzo and R. Gray, "An Algorithm for Vector Quantizer Design," in IEEE Transactions on Communications, vol. 28, no. 1, pp. 84-95, Jan 1980.

It was again extended in 1989 for entropy coded vector quantization:

P. A. Chou, T. Lookabaugh, R. M. Gray, “Entropy constrained vector quantization,” IEEE Transactions on Signal Processing, vol. 37, no. 1, pp. 3142, January 1989.

I can recommend visiting www.data-compression.com/vq for a descent explanation of the scheme. It works in a simple way:

Divide the vector space.
Update the partitions.
Repeat the procedure for the new partitions until the quantization error is acceptable.

At the end, all data points that belong to a partition are quantized to the coordination of the center of partitions. The centers are indexed such that to represent all quantized values only the index is stored/sent.

Animated illustration from the same website for 2D Gaussian source:


How to specify the input data bit width for fixed point implementation in windowed OFDM?

Depends on which processor platform you want to implement it. If it is an ASIC (Application Specific Integrated Circuit) chip, they come with ARM processor cores in most cases and real time computation require fixed-point implementation. There is also memory limitation in such processors due to smaller footprint of device. If you implementation is on processor which is 32 or 64-bit wide and has large amount of memory at its disposal (such as on Laptop, Desktop), you can use 32-bit fixed point notation (even though floating point would just work fine like IPP DSP Library from Intel). In your case, if your I and Q data each are 16-bit fixed point numbers, I don't think you have to worry too much. I am not sure if there is a specific advantage but almost all fixed point convention I have come across use 1 sign bit and rest fractional bit. In this case Q1.15 format. In that case you need to normalize all your symbols, channel co-efficients etc so that they are within +/-1.


Why does the periodic signal in time always give a discrete frequency spectrum

Here's an intuitive explanation if the convolution theorem is taken for granted:

Since the time-domain signal is periodic, one can say that it can be built by "copying and pasting" the same block of signal every period: your periodic signal can be expressed as a little block of signal (spanning one period) convolved with a dirac comb.

Thus, its Fourier transform will be the Fourier transform of the little block multiplied by the Fourier transform of a Dirac comb (which is another Dirac comb). Multiplying a continuous signal by a Dirac comb yields a discrete signal.

The same reasoning is also true the other way round (discrete in time implies periodic in frequency).


What can we deduce about variance when we are given the noise spectral density?

Given a low-pass signal $x$ that is a realization of random variable $X$ with constant spectral density $a=N_0/2$ in $f=-B\ldots B$ (and zero otherwise) we can calculate the mean power of that signal by $\bar P = 2aB=N_0B$, where $B$ is the signal bandwidth. The mean power of a random variable is defined as the expectation of $X^2$, i.e. $\bar P = \mathrm E[X^2]$. The correspondance between mean power and variance $\sigma_X^2$ of a random variable $X$ is
$$
\sigma_X^2=\mathrm E[X^2]-m_X^2
$$
where $m_X$ is the mean value of $X$. That is, for zero mean (which is often true for noise) we get
$$
\sigma_X^2=\mathrm E[X^2] = N_0B
$$

Concerning your second question, I assume that $T_0$ is the symbol duration. Thus the symbol rate $R_\mathrm{S} = 1/T_0$. The bandwidth of a signal with symbol rate $R_\mathrm{S}$ is often approximated by $B=R_\mathrm{S}/2$ and thus
$$
\sigma_X^2 = N_0B = N_0R_\mathrm{S}/2 = N_0/(2T_0)
$$
(You can also check the units in your equation: $N_0$ is power per Hz and $T_0$ is in sec. But the result should be a power)


How Is Laplacian of Gaussian (LoG) Implemented as Four 1D Convolutions?

I think Chris Luengo's answer is perfect.
The trick is that you can calculate the 2nd derivative of the image (Using Finite Differences -> Convolution) and then blur it with Gaussian Filter.
Since both are seperable kernels you can do that by 4 1D convolutions.

If the input image is given by I

Calculate Ixx (The 2nd derivative on x direction) using convolution.
Filter Ixx with 1D Gaussian Kernel along the x direction.

Do the above for the y direction as well.
The LoG image is the sum of both.

To answer your questions:

See my explanation above. I think it matches your understanding.
This is the very definition (The Discrete) of the LoG.
I don't see the contradiction between #62 and #63.


Digital beamforming and phase shift stage insertion

In a first analysis, it does not matter where the phase shifting is inserted, as long as the mixing stage is linear.

Indeed, shifting a temporal signal is equivalent to multiplying it with a complex number integration the phase shift $\theta$ :

Unshifted signal : $s(t)$
$\theta$ phase shifted signal : $s'(t) = s(t).e^{j\theta}$
$f_0$ carrier mixed unshifted signal : $s_c(t) = s(t).c(t) = s(t).e^{2j\pi f_0.t}$
$f_0$ carrier mixed $\theta$ shifted signal : $s'_c(t) =s(t).e^{2j\pi f_0.t}.e^{j\theta}$ = $s(t).c'(t)$

where $c'(t) = e^{j(2\pi f_0.t+\theta)}$ is the shifted carrier.

In the end phase shifting the signal to be modulated is equivalent to phase shifting after carrier modulation.


Deep Learning: Classification vs. Convolution for Signal Restoration

By deep learning, I'll assume you mean neural network. To develop a neural network, you'll need labeled data. This means you need a bunch of example inputs ($r$) and outputs ($X$). Once you have that, you need to build a network. As far as I can tell, currently the guidelines for building these things are very loose and ad-hoc. There is general consensus on good or bad activation functions and general suggestions, but there is no formal theory about how many neurons to choose, how many layers, etc. (at least since the last time I took a dive into this).

Common problems will arise such as overfitting and there are techniques to help this but you need to be smart enough to spot the signs and then apply the fixes. For example, you can read about how observing your training and testing accuracy can help determine if a network is overfitting.

The better part of this answer is this: what you're describing is called deconvolution (https://en.wikipedia.org/wiki/Blind_deconvolution) and has been accomplished using "traditional" signal processing techniques long before the deep learning craze. Personally, any time a problem has a structured solution, I go in that direction. Reliable neural networks can be hard to develop despite all the helpful libraries in the world and when it comes to troubleshooting I think you'll find more traditional "tried-and-true" methods to be much easier to not only get working but also to understand why they work and this will lend itself much better to be able to adapt to your specific problem.


I prefer to think of it as expanding the spectrum due to increasing the sample rate.  When you expand the frequency spectrum you reveal the aliases that were already there in the other Nyquist zones.  See this answer for more detail about this.


Is there a rule of thumb for selecting the variance of the input?

You can define your model as $$\hat{x}[n] = \frac{1}{||\mathbf{x}_{n}||}\mathbf{A}^T\mathbf{x}_{n} +\frac{1}{\sqrt{\mathsf{SNR}}} w[n]$$
In this way, you define a range for $\mathsf{SNR}$ in dB (e.g. between 0 to 10 dB), then translate it to natural (non-dB) scale, and then $\frac{1}{\sqrt{\mathsf{SNR}}}$ gives you the noise power. The signal power is normalized to one which makes it easy.

Here, $w\sim \mathcal{N}(0,1)$ is a standard Gaussian RV and you can simulate it as w = randn;. In fact, $\frac{1}{\sqrt{\mathsf{SNR}}}w=w'\sim \mathcal{N}(0,\sigma^2_{w'})$, where $\sigma^2_{w'}$ has the right proportion wrt the SNR.


Failed to Get Silence in Audio Streaming

the PCM data consist of signed value between 0-255.

That seems to be your problem. A 16 bit MP3 file decodes into signed integer values of -32768 to +32767. Perhaps you are looking at bytes and not at samples ?


How to verify if two images have the same object/content?

I am no expert, but I would like to give some pointers which could help you to solve.

Regarding the method you have mentioned/used, is for image registering/stitching. You are interpreting the results wrong. Calculating the percentage of feature match is not the intention of inliers, rather computing corresponding feature points in both images which could be used for image registering/stitching. Inlier is the feature point for which the algorithm has found a corresponding feature point in the other image. Feature points for which the algorithm has not founda match are called outliers. More inliers, better the match. Keep in mind, the algorithm can sometimes match some feature points of completely unrelated images(try running your method on two unrelated images)
Can you use the method you have mentioned to solve your problem ?.
Yes,You can use number of inliers for each image in the database and decide the one that has the highest number of inliers as the match. You can have a minimum count of inliers to be present, in case the picture of a video-game doesnt exist in your database.
Other ways/methods which you could explore to solve this.

Explore more on the features you can use. SURF does well for some applications, for some it doesnt work well. You have to choose based on your data (image database). You can check with SIFT, and also color histogram as features.
You can check bag of words for image search. (check this)
Whenever you add a new game to database, some algorithm(eg - bag of words) need the code book to be recomputed everytime you add a new image to database, so choose the algorithm wisely , coz as the image database gets bigger, the cost of computing  increases.
Also think in lines where scale of the images differ(i.e image is database and search query image) Some features like SURF/SIFT are scale invariant, if you choose to use any other, keep this in mind.

These points might help you in deciding and I hope I did not complicate your life :)


FFT and IFFT are linear operators, and as such, the results only make a lot of sense in a linear intensity space, not if indexed into a non-linearly mapped space.


Non-symmetrical lemniscate curve parameterization

There is a polar asymmetric modelling for the lemniscate:

which you can rotate and scale more easily that in Cartesian coordinates (but easy to convert). Such parameterizations are used very often as they can be less troublesome to fit.

Similar curves in astronomy are also called analemmas:

.

From an image processing fit point of view, you can consult: A Unified Scheme for Detecting Fundamental Curves in Binary Edge Images, Asano, Tetsuo et al., Computational Geometry, 2001


Removing transients in highpass filtering with MATLAB

Your first order filter recursion for some real constants $a,b,c$ is

$$ y[n] = a x[n] + b x[n-1] - c y[n-1] $$

with the two initial memory states $x[-1]$ and $y[-1]$ at $n=0$.

Your "no transient" condition can be translated to $y[0]=0$ and a necessary second condition so that you can solve for both of your memory states. The second condition could be, that the discrete derivative of $y$ also vanishes at $n=0$, so $y[0]-y[-1]=0$. You can also take any other condition that seems sensible to you.

The two equations give you a unique solution for the two unknown memories, namely:
$$y[-1] = 0$$
 and
$$x[-1] = -\frac{a}{b}x[0] $$

Alternatively, your conditions may better be chosen as $y[0]=0$ and $y[0]-y[-1]=x[0]-x[-1]$ in order to capture the initial slope of the input. The resulting recursion equation at $n=0$ is then
$$0=a x[0]+b x[-1]+c(x[0]-x[-1])$$
giving you the solution
$$x[-1]=-\frac{a+c}{b-c}x[0]$$
and
$$ y[-1]= -\left(1+\frac{a+c}{b-c}\right)x[0]$$

(Please check my calculations!)

But in general you cannot expect a simple initial condition to give you the same result as knowing the signal history. So you can only take this to a certain point and in general it would probably be better if you just discarded the transient response of your output.


Just as the DFT assumes a periodic signal by construction, the DCT assumes an even signal by construction. We can think of taking a DFT of a non-periodic signal by extending it periodically, but restricting our "viewing window" to the original signal itself. Now, what happens when you take the DFT of an even signal? The DFT is then real (which means it can be constructed with only cosines as basis functions). So, in a similar spirit to the periodic extension, we can take an arbitrary non-periodic, non-even signal, mirror it, then extend it periodically. This results in a real, cosine-based DFT. Again, we restrict our "viewing window" to the original signal itself. Note that although we've doubled the length in the frequency domain, we only have real components. Effectively this means we store the same amount of information compared with the traditional DFT. The reason the DCT is useful in compression is that compared with the traditional DFT, there is more energy in the lower frequency components than in the higher frequency components (energy compaction), so we can compress those higher frequency components. This is true because, as you said, this extended mirrored periodic signal lacks the sharp discontinuities of the traditional DFT input and as a result is smoother.


how to calculate the maximum obtainable directivity for a microphone array?

When using a microphone array for baseband acoustic signals, Array Gain is a preferred performance metric over Directivity Index.

In geophysics, the term “stacking” is often preferred over the equivalent “delay-and-sum beam former. It is certainly more descriptive because you essentially time align signals like a stack and add them coherently.

If the noise is independent at each microphone, for the aligned stack, the upper bound on coherent gain is $10 Log N$.

Directivity is a function of the beam pattern which is a function of a single frequency.  For an RF system with an obvious carrier, the pattern at that nominal frequency is an often used way to characterize gain.

The stacked array gain is baseband. It is the upper bound on gain.


audio signal time delay estimation(TDE) in indiviual and separate mics

single mic, i'm not sure.  it seems to me that all you can do is mark different events (perhaps with envelope following or perhaps by comparing the signal or parameters measured from the signal to templates) and just use the time difference of the onset times of each marked event.

for two mics, there is cross-correlation.  but some people have also suggested an adaptive LMS filter (one of the mics is the "desired signal", $d[n]$) of some sort and computing some measure of group delay from the FIR coefficients that result in the LMS.  something in me makes me think is the Julius Smith/Perry Cook thing, but it might be someone else.


Why only odd harmonics after non-linear amplification?

An ideal amplifier would have a transfer characteristic of $f(x)=Ax$: the input signal comes out amplified and otherwise undistorted. A real amplifier will deviate from this and go into saturation. We could model it by a polynomial $f(x) = \sum_{n} a_n x^n$. Now, what we would still expect is that the amplifier treats positive and negative values the same way, i.e., it is symmetric in the sense that $f(-x) = f(x)$. This leads to the requirement $a_n = 0$ for all even $n$.

Now, one can show that raising a sinusoidal with frequency $f_0$ to the $n$-th power generates harmonics up to order $n$, however, preserving their oddity:

$n = 2$ generates frequencies 0 and $2f_0$
$n = 3$ generates frequencies $f_0$ and $3f_0$
$n = 4$ generates frequencies 0, $2f_0$ and $4f_0$
$n = 5$ generates frequencies $f_0$, $3f_0$ and $5f_0$

and so on. It's a consequence of the binomial theorem and quite easy to show (see, e.g., here).

In short: A transfer function that is odd symmetric will generate only odd harmonics. This is why you are seeing only those.


Coding a Variable Delay with LFO (MATLAB/C)

I answer this way because i don't have permission for comments.
I don't have a complete solution. What happens with an chorus flanger is this.

Signal xt goes in split in two(or maybe more) signal paths and mixed together at the end.
The difference is the length of the delay (Flanger <10ms and chorus 20 - 60ms).
First your Lfo discription isn't right.
It should me something like :lfo = A *sin(2*pi*lfo_rate*[0:length(output)]/fs) + Offset; Where A is half the time (Positive and negative) and offset the center time where the LFO is turning around.
And `A =< offset' otherwise the time gets negative

The difficulty is i think a variable buffersize. let say a chorus center time is 50 ms with ±10 ms delay. no Lfo you have 44100 * 20/1000 = 2205 samples. the max buffer = 2646 samples and the min buffer = 1764 Samples.

Hope you can do something with the info


Equalizer coefficients and channel coefficients

With a blind equalization technique like the constant modulus algorithm (which is often implemented using a least mean squares (LMS) filter as you indicated), you aren't directly estimating the channel impulse response itself. Instead, the signal model is like this:

The receiver observes the following signal:
$$
x[k] = s[k] * c[k] + n[k]
$$

where:

$s[k]$ is the underlying desired signal,
$c[k]$ is the channel's impulse response, and
$n[k]$ is a noise process introduced by the channel (typically assumed to be additive white Gaussian noise)

A blind equalizer operates by taking advantage of some constraint on the properties of the desired signal $s[k]$. For instance, if it is constant envelope, like phase-shift keying, then the constraint could be:
$$
|s[k]| = 1
$$
The CMA-LMS blind equalizer operates in two steps. First, it uses the current length-$L$ vector of LMS filter taps (where $L$ is a design parameter of the filter) $\mathbf{w}[k]$ to calculate its output:
$$
y[k] = \mathbf{w}[k]^H\mathbf{x}[k]
$$

where $\mathbf{x}[k]$ is a vector of the last $L$ input samples to the filter. Next, it uses the calculated output to update its set of taps for the next iteration:
$$
\mathbf{w}[k+1] = \mathbf{w}[k] - \mu \mathbf{x}[k](|y[k]|^2 - 1)y^*[k]
$$

When compared to the typical LMS filter structure, one can identify the error signal $e[k]$ as:

$$
e[k] = (|y[k]|^2 - 1)y^*[k]
$$

You can see here how the error captures how far the filter output signal $y[k]$ deviates from the expected constant envelope of the desired signal $s[k]$.

Looking at this structure, you can see that at no point is the equalizer trying to directly estimate the channel's impulse response. Instead, it's more "given the knowledge that I have of the expected signal's structure, how do I need to adjust my filter taps to enhance that structure in the filter output?" So it's more of an inverse problem; the tap vector will tend to more closely resemble the inverse of the channel impulse response. Blind equalization techniques like this can do a surprisingly good job at doing so (even more surprisingly, techniques like CMA can even work for signal types that aren't constant envelope, like QAM).


Effect of changing sample rate, window duration and zero padding on DTFT and DFT

halving the $F_s$ (keeping $T$ and $N$ constant) would remove the right half of the bins (halving the maximum frequency)

No. Number of bins remain same because you have not changed $N$. What will change is the resolution. Each bin will represent $F_s/(2N)$ instead of $F_s/N$. Your DFT will not represent analog frequencies from $F_s/2$ to $F_s$ That is why you misunderstood as 'remove right half of the bins'.

  halving $T$ would change the shape of the DTFT, halving the number of lobes. Without zero padding this would halve $N$, but $F_{max}$ would remain constant and so every second bin would be removed

I am assuming that when you halve $T$, you are also reducing DFT to $N/2$. If you reduce FFT size with no zero padding, yes, you are dropping every second bin. Because the resolution now increases to $F_s/(N/2)$. For a period signal, with $N/2$ multiple of period, you will still see only 2 peaks for the FFT because the lobe width also has doubled, and the zero crossing appear exactly at multiples of $F_s/(N/2)=2 F_s/N$.

  doubling N by zero padding the DFT would keep the DTFT and Fmax constant but would double the number of bins

Yes, zero padding would increase the resolution of DFT. $F_s/(2N)$ is the resolution as compared to earlier $F_s/N$. But remember by zero padding you are not adding any new information about the signal. It just reveals more of the DFT you already have.


Unfortunately the only way I know of to convert their .ptw format to an standard video format is by using FLIR's ResearchIR software:  http://www.flir.com/thermography/americas/us/view/?id=51371

ResearchIR also supports exporting to .mat files too.

Hope that helps,
Andy


N-connectedness means that from any pixel to any other within the component, there should exist a path composed of N-connected steps. In other words, as long as two groups have one connected neighbor, they should be merged into a single group. It does not mean a 2 or 3 pixel group, but rather a component which consists of linked pixels.

See the connectedness garph theory.

With such definition, it is obvious that the author is right.


Convert RGB Image to Grayscale and Display It (Python + Matplotlib)

You need to set the cmap properties of the imshow() function:

https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.imshow.html

Use:

plt.imshow(img, cmap = "gray")


Confusion in Peak to Average Power Ratio (PAPR)

Preference for low PAPR is subjective. From a power efficiency and unit cost perspective, low PAPR is certainly beneficial, but from a spectral efficiency perspective, it typically isn't. High entropy modulation means high PAPR (i.e. signal looks more like Gaussian random process).

A low PAPR can mean that less backoff is needed in a power amplifier, and hence better power efficiency can be achieved. Similarly, low PAPR means less headroom required in an ADC and thus allows for lower-precision in the ADC. These back-off/headroom requirements for TX and RX are closely related to PAPR. Some typical systems considerations that are linked to system requirements:

modulation - higher-order QAM and/or multicarrier modulation means higher PAPR, which means an increased budget for headroom across signal paths.
fading - not normally a part of PAPR, backoff for fading is closely related and arguably becomes an overall PAPR consideration for RX paths (we normally reference PAPR as signal-related, but fading typically increases the ratio)
code rate and blocksize - high coderate and large blocksize, particularly at high-order modulation, means that the signal can be highly susceptible to clipping, non-linearity, and other impairments, and thus we need to allow headroom for the rarest of peaks (definition of PAPR is somewhat subjective, too; are we worried about the 1E-3 probability peak or the 1E-6 probability peak?)
bitwidth - ADC, DAC, and processing chain bitwidth. Higher PAPR plays here, too. For higher PAPR, more backoff from full scale is needed to avoid distortion of the peaks. This requires more bits and analog linearity for signal representation and processing.

For multicarrier modulation (OFDM, DMT, etc), PAPR can certainly be problematic. To visualize, try modulating a high-order QAM constellation on subcarriers using an FFT (IFFT). The sum of tones results in a constructive and destructive combining across the IFFT symbol, giving high PAPR.


Understanding the mathematical proof for the alias frequencies in a sampled sine wave

The reason is that if it is true for any $m$, it is also true for $m=kn$.

I will sketch the proof in another way.

Call $f_s = 1/t_s$ sampling frequency where $t_s$ is sampling period, the two signal $x(t) = \sin(2\pi f_0 t)$ and $x_k(t) = \sin(2\pi (f_0 + k f_s) t)$ have the same values at sampling instants (aliasing), i.e. $x[n] = x_k[n]$.

Indeed,

\begin{align}
x[n] &= \sin (2 \pi f_0 n t_s) \\
x_k[n] &= \sin \left(2 \pi (f_0 + k f_s) n t_s\right) \\
&= \sin (2 \pi f_0 n t_s + 2\pi k f_s n t_s) \\
&= \sin (2 \pi f_0 n t_s + 2\pi k n) = \sin (2 \pi f_0 n t_s)\\
&=x[n]
\end{align}


How the noise in time series can be reduced?

Yes indeed, there is a huge body of works on the topic of variational methods for signal restoration (more general than denoising), one example being total variation denoising, or proximal methods.

If you can be more specific, you could get more accurate and precise answers.


The first thing to note is that the LMA algorithm is an estimator, and tries to minimize error in the least squares sense, between your data, and a model function of variable parameters.

Typical of LSE based frequency estimators, you need to seed them at values close to the global minima of the cost function, otherwise they will find another local minimum to lock on to, and hence they need to be seeded properly, or 'in the ball park', ostensibly using some apriori information. How to attain this a-priori information? See below.

  This solution worked wonderfully for toy data that I generated, but it
  fails quite frequently on real data, because of the noise. The
  amplitude of the noise can be as large (at most) as the amplitude of
  the signal.

LMA will provide you with an amplitude, phase, and frequency estimate in the LS sense. Of course, the more noise there is in your data, the worse off this estimate is. This is why you will need to pre-process your data to remove as much noise as possible, before the application of the LMA.

You mention zero-padding your data, and while this will not increase your frequency resolution, this will certainly increase your frequency granularity. (Interpolation in the frequency domain). With those facts, I would attempt the following solution:

Pre-Process the data via an (FIR) Band-Pass Filter centered around my frequency of interest.
Perfrom a zero-padded FFT so as to ascertain the peak frequency of interest within the BPF'd data's band.
From this peak, ascertain the frequency, amplitude, and phase, with which to seed LMA with. (Note, make sure to remove the phase offset the BPF filter added for the seeding value).
Run the LMA.

  I'm guessing you mean throw away the highest frequency components etc?
  Would the statistics of my noise matter in that case? I can assume
  it's gaussian, but that may not be strictly true in all cases

The statistics of the noise does not matter in this case. All you are doing is removing frequency bands where noise exists, independent of their statistics. In this way, you are not utilizing DFT bands where your signal does not exist, and retaining the DFT bands where your signal exists.


Odd problem with really long sinc filters

A frequency sweep has non-zero bandwidth over the length of any FIR filter.  Use a long enough filter (compared to the sweep rate) and nearly half the spectrum of the sweep  could end up on each side of the cutoff.

Also, use a narrow enough transition band, and you might end up with Gibbs horn ringing right near the edge of the transition frequency.


What's wrong with this Whittaker-Shannon-Kotel’nikov interpolation implementation?

Replace the line

y+=samples(k)*normalized_sinc((t-k*T_s)/T_s);

by

y+=samples(k)*normalized_sinc((t-(k-1)*T_s)/T_s);

because the indices inside the sinc function must start from $0$, not from $1$. The first sinc must have its peak at $t=0$, not at $t=T_s$.


What are the units of amplitude of my ultrasonic recording equipment

Given that the gain of the system and/or ADC seems to be unknown, the scale of the sample values can only be determined from samples of a calibrated source of known value at a known distance from the transducer.  With a variable enough gain knob (or AGC) the scale of the samples could be anything.


As pointed out in the comments, this question is rather broad, so here I will just answer your last question about the autocorrelation of the output signal $y[n]$. You might know the following relation concerning the power spectra of input and output sequences:

$$S_y(\omega)=S_x(\omega)|H(\omega)|^2\tag{1}$$

where $H(\omega)$ is the frequency response of the LTI system. From the inverse Fourier transform of (1) you get a corresponding relation between the autocorrelation functions:

$$\phi_y[m]=\phi_x[m]\star h[m]\star h^*[-m]\tag{2}$$

where $h[m]$ is the system's impulse response, and $\star$ denotes convolution.


Re-mixing high and low pass filtered audio has interference in the shared transition band?

You are assuming that a parallel high pass and low pass add up to "flat". In most cases, they do not. You have to carefully design both the phase and the amplitude responses in the transition band to match each other.

There are a few classes of filters that are useful here

1st order lowpass/highpass to add up to unity
Odd order butterworth filters add up to an allpass filter (flat frequency response, but some phase distortion). The two paths are either 90 or 270 degrees phase shifted (depending on order)
Linkwitz Riley filter (two butterworth cascaded) also add up to an allpass (when choosing the appropriate sign for recombination). The two paths are in phase or out of phase (depending on order).
linear phase FIR filters that are designed complementary.
Specific types of filter banks (reconstructing)

All of this only helps if your processing doesn't induce significant differential phase shifts. Unless this phase shift is constant and can be compensated, you will always get some amount of interference.