Detecting fundamental zero crossing point I am looking for ideas from the more experienced members on algorithms for determining the fundamental zero crossing point of a 50 or 60 Hz signal. The fundamental zero crossing point is then used as the commencement point to calculate the RMS value of the original (voltage) signal AND other (current) signals.

Typically the frequency content of the Voltage is below 5% THD. At the moment I can't give u an exact figure on the SINAD but the typical Voltage signal will almost always be at 75% full scale (~250mV), or close to, of a delta sigma ADC (SINAD 110db) channel, with a sampling rate set to 6.4kHz. The target is an embedded DSP ARM M7 ~200Mhz with floating point hardware.

Ideally i would like to have 3-4 algorithms to prototype in Matlab and find the most appropriate. FIR vs. IIR vs. Bandpass vs. Zero-phase? Smarter approach? Would like to understand issues relating to sensitivity to noise, harmonics, group delay, edge effect, etc...

The IEC 61000-4-30:2015 standard on Power Quality states:
"r.m.s. voltage, refreshed each half-cycle
value of the r.m.s. voltage measured over 1 cycle, commencing at a fundamental zero crossing, and refreshed each half-cycle.

r.m.s. current, refreshed each half-cycle
value of the r.m.s. current measured over 1 cycle, commencing at a fundamental zero crossing on an associated voltage channel, and refreshed each half-cycle."

Thanks
Alex

I am not confident I have the best algorithm, but short of other suggestions I could offer you this idea:

Use an NCO (Numerically Controlled Oscillator) clocked to a known and accurate time reference to create a 50 Hz or 60 Hz reference signal, and lock the NCO to the power line signal using a Phase-Locked-Loop (PLL). The NCO will lock to the power line signal in quadrature phase with optimum noise suppression (based on the loop design and desired time averaging versus dynamic requirements). By reading the phase accumulator word the instantaneous phase offset between the power line signal and the time reference can be determined.

This idea is developed as follows:

If you have a known time reference and can generate a reference 50 Hz or 60 Hz sinusoidal signal, then you can use a phase detector (by multiplying the two signals and low pass filtering) to determine the phase difference between the two, and from the determined phase can determine the zero crossing relative to the known time reference. There will be calibration involved to offset the time delay of the low pass filtering. The advantage of this approach over single sample decision metrics (threshold zero crossing detectors) is its relatively high noise immunity based on the frequency cut-off (averaging time) of the low pass filter used.   The phase detector is based on the mathematical relationship between the product of two sinusoids at the same frequency; the result is an output magnitude that is proportional to the cosine of the phase of the two signals. This provides an unambiguous mapping of phase to magnitude over a 0° to 180° range.  Most significant is the result is averaged (via the low pass filter) over all samples taken, so benefits from significant noise averaging- EVERY sample in the waveform is used and contributes to the estimate of the waveform's phase position!

This is not yet complete as the primary drawback with this approach is that if the two signals approach 0° and 180° phase, the slope of the phase detector output approaches 0; meaning the result is not very sensitive (accurate) to phase shifts. The ideal operating point is when the two signals are in quadrature as the output would have the highest slope (sensitivity) to phase changes.
This leads to the improved approach of using a reference signal where you can precisely control the phase, and adjust the phase to the quadrature relationship between the signals by monitoring when the output of the phase detector is zero- this is basically how many phase locked loops are designed to operate, and you could use such a phase locked loop (PLL) approach to do this adjustment automatically; as long as you have accurate knowledge of the relationship between the control to your reference waveform and its phase. Further this gives you well designed control of your averaging time via the loop bandwidth of the PLL design. A numerically controlled oscillator is such a device where you have precise knowledge of its phase at any given time, based on its phase accumulator.

For more details on the NCO implementation, see this link: Numerically Controlled Oscillator (NCO) for phasor implementation?


Scale signal to a given range Is there a function in MatLab to scale a given signal to a specific range, like how soundsc does it with the input?

I have to signals, one is a wave file, using wavreader (I know one should use audioreader now) to load the file, and another file that is being read using fopen and fread.
Now I want to compare those to signals, but the values of wavreader seems to be between -1 and 0, but the other file's signal goes anywhere between -5000 and 0.

Given the signal $\mathbf{x}$, you can scale the values of $\mathbf{x}$ to fall on the interval $[a,b]$ with

x_scaled = (x-min(x))*(b-a)/(max(x)-min(x)) + a;

Be aware that this is a very naive scaling, but it is the best you can do without any extra information.  When I say naive, I mean that it assumes that the scaled signal should always have a minimum value of $a$ and a maximum value of $b$.


"Flanging" in frequency domain? Flanging is defined as a mix of two identical signals where one signal is delayed in time by a small and gradually changing period, around 10-20 milliseconds.

Since delay in the time domain is equivalent to a phase shift in the frequency domain, could then a Flanging sound effect be done in frequency domain?

Could we do a DFT transform, make a copy of amplitudes and phases, shift phases (and gradually change the shift between -180 and +180 degrees), and add the result to original amplitudes and phases?

Frequency domain does not really work here. A flanger is NOT a linear time invariant system. You could try to implement this is "piece-wise" LTI and change the delay only at a frame boundary but this will mostly likely result in unacceptable artifacts.

A flanger will require a time variant fractional delay. Rounding the delay the nearest integer delay results in unacceptable "zipper" noise.


What are the limitations of FIR filters? I know that for type II FIR filter it adds a zero at $z=-1$ so it cannot design a high-pass filter.

My question:

  what type of limitations the other types have (type I, III, IV)?

Note that this is not about general limitations of FIR filters, but about the special case of linear-phase FIR filters. If you understand why a type-II linear-phase FIR filter has a zero at $z=-1$, then the limitations of the other types should be obvious too.

It's always about zeros at either $z=1$ (DC) or $z=-1$ (Nyquist). Given the transfer function

$$H(z)=\sum_nh[n]z^{-n}\tag{1}$$

it is easy to see that

$$H(1)=\sum_nh[n]\tag{2}$$

and

$$H(-1)=\sum_n(-1)^nh[n]\tag{3}$$

Now look at examples of the four linear-phase FIR filter types. An example of an impulse response of a type-I filter is

h = [1 -1 2 3 2 -1 1]

Using $(2)$ and $(3)$ can you say anything about zeros at $z=1$ or $z=-1$?

An example of a type-II filter is

h = [1 -1 2 2 -1 1]

Clearly, applying Eq. $(3)$ results in a zero at $z=-1$, regardless of the actual values of the impulse response.

The following two filters are type-III and type-IV filters, respectively:

h = [1 -1 2 0 -2 1 -1]

h = [1 -1 2 -2 1 -1]

Now apply Eqs $(2)$ and $(3)$ and see what you get.

This answer contains more information on the four types of linear-phase FIR filters.


How to cascade two filters? I have to filter an audio signal through a cascade of two equal low-pass FIR filters. They are half-band filters so each one has a cut-off at $\pi/2$.

Then, I have to multiply the resulting signal (obtained through the cascade of the two filters) for $z^{-q}$, where $q$ is a scalar.

Firstly, I thought to convolve the impulse responses of the two low-pass filters in this way  (where hlp1=hlp2):

hconv=conv(hlp1, hlp2),

and then filter the audio signal (called x) in this way:

x_filtered=filter(hconv,1,x)

Is it correct?

I expect that the frequency response hconv is a low-pass filter with cutoff at $\pi/4$, but the command fvtool(hconv) shows that hconv has a cutoff at $\pi/2$, and not at $\pi/4$.

Is it correct?

Moreover, how can I implement the multiplication of the filtered signal x_filtered for $z^{-q}$?

You can easily verify yourself that the following two ways of cascading two filters are equivalent (up to numerical inaccuracies):

% method 1
y1 = filter(h1,1,x);
y = filter(h2,1,y1);

% method 2
h = conv(h1,h2);
y = filter(h,1,x);

So you were on the right track when you convolved the two impulse responses. Note that linear convolution is what you need, not circular/cyclic convolution.

As for the multiplication with $z^{-q}$, where I hope that $q$ is an integer, you need to understand that it is not the (time-domain) signal that is being multiplied, but that this multiplication refers to the Z-transform of the signal. In the Z-transform domain, a multiplication with $z^{-q}$ corresponds to a delay of $q$ samples.

Finally, regarding your expectation that the cut-off frequency of the cascaded system is $\pi/4$, note that the frequency response of the cascaded system is simply the multiplication of the two individual frequency responses, so the cut-off frequency does not change significantly (it may change a bit depending on how you define it). The transition from pass band to stop band will be steeper, and the stop band attenuation increases, but the start of the transition band remains unchanged. If both filters are identical you will basically get a better low pass filter with the same pass band and stop band edges.


How does one calculate a pole-zero plot? To my understanding, pole-zero plots are used to analyze or visualize transfer functions. Suppose there is some very simple system, for example a simple low-pass filter (so it is linear and time-invariant). How does one calculate the pole-zero plot of such system?

Suppose you are given a system with transfer function

$$H(z)=\frac{(1-3z^{-1})(1-7z^{-1})}{(1-4z^{-1})(1-6z^{-1})} $$

Poles

Poles are the values of $z$ for which the entire function will be infinity or undefined. So, they will be the roots of the denominators, right?

Look here, what values of $z$ will turn the transfer function tend to infinity? Obviously it's $z= 4$ and $z=6$, because if you let $z$ equal 4 or 6, the denominator will be zero, which means the transfer function will tend to infinity.

Zeros

Zeros are the values of z for which the transfer function will be zero. As you have guessed correctly, zeros come from numerator. In this case, zeros are $z= 3$ and $z=7$, cause if you put $z= 3$ or $z=7$, the numerator will be zero, that means the whole transfer function will be zero.

So here poles are $z=4$ and $z=6$, and zeros are $z=3$ and $z=7$.

Pole-Zero Plot

About finding the Pole zero plot, you draw a complex plane. Then you put the values of poles as 'X' marks and zeros as 'O' marks. Here I took the liberty of drawing the pole zero plot of the system:

So, for low pass filter, you find out the transfer function, then the poles and zeros.


Proving linearity and time-invariance This question was asked in one of our quiz. Question is

Determine if the system with the input-output relationship given by
$$y(t) = \int_{t-10}^{t} \cos(\tau)x(\tau)d\tau$$
is linear and time invariant.

I was able to prove linearity but was not able to prove the time invariance of the system.

My Try
For invariance we need to show output resulted from $x(t-t_o)$ is same as $y(t-t_o)$.
Let $$y_1(t) = \int^{t}_{t-10} \cos(\tau)x(\tau-t_o)d\tau.$$
Now when I change of variable to $z=\tau - t_o$, $dz = d\tau$ it leads to
$$y_1(t) = \int^{t-t_o}_{t-10-t_o} \cos(z+t_o)x(z)dz$$ where I get stuck (using $\cos(A+B)$ identity doesn't seem to work well) and not able to proceed further.
For $y(t-t_o)$ integral becomes $$y_2(t-t_o) = \int^{t-t_o}_{t-10-t_o} \cos(\tau)x(\tau)d\tau$$
From here it seems to time variant system. But when $t_o=2\pi n$ it becomes time invariant. But system cannot be time invariant and time variant. That's the doubt I have with this question whether the system is time invariant or time variant and why.

Your proof is correct, and the result is that the system is time varying because the response to $x(t-t_0)$ does not in general equal a delayed response to the input $x(t)$. Of course, for the special case $t_0=2\pi k$ the delayed response to $x(t)$ equals the response to $x(t-t_0)$, but for time-invariance this equality must hold for any $t_0$.


autoregressive moving average code implementation I am new to DSP and i am trying to take a wav (human speech) file and apply ARMA filter and plot its PSD graph in python.

I see that there are a lot of AR implementations but almost none ARMA. I found only statsmodel that implemented ARMA but it requires the data to be stationary (but sound with lags of 30 is almost stationary)

Does anyone have a working piece of code that i can use to plot?
If not, i will be happy to get some reference on how to implement it.

Thanks in advance!

I dont know in Python, but this simple matlab/octave loop code:

function y=arma(x,A,B)
na=length(A);
n=length(x);
y(1:na,1)=zeros(na,1);
for i=na:n
    y(i,1)=1/A(1)*A(2:na)*y(i-1:-1:i-na+1)+1/A(1)*B*x(i:-1:i-na+1);
end

Implements a ARMA filter, with polynomials $A$, $B$, such as per this 5th order filter:

A=[1 0 0 0 0 0];
B=[0 0.2 0.2 0.2 0.2 0.2];
n=1000; dt=0.001;
t=dt*(0:n-1)';
x=randn(n,1);
y=arma(x,A,B);
plot(t,[x y]);


Introduction to Clustering and Classification of Linear Sequence I'm a programmer/engineer that is getting into developing prototype wearable technology. I would like to be in a position to interpret raw accelerometer data and categorised the data into activities like running/walking/rest, etc.

The building and coding of the hardware and software is relatively straight forward, but I do not know the terminology and mathematical processes required in order to take the raw data and turn it into a training set that can be used for a real time system.

My question is this: What is a good DSP starting point to take a time based linear series and extract features so that they may be identified at a later stage?

At present, I have an tri-axis accelerometer device which samples at 50Hz with a range of -/+4g. I would like the capability of being able to readily identify activities based upon a training set, or machine learning process, if that's even possible?

Apologies that this isn't necessarily a specific DSP question, thanks for you assistance.

This seems very much like a fit for Hidden Markov Models (HMM). They are pretty powerful in handling sequence based input such as time series.


Reducing standing waves with phase randomization Will constantly randomizing the phase of an audio signal, say a speaker in the corner of a square room, reduce standing waves (i.e. room modes) in the room?

For example if you wanted to create a diffuse field (i.e. no standing wave interference patterns) in a small room at low frequencies.

Assuming that you are not changing the phase too much, you are basically turning sinusoids to narrow-band noise. If there is a room frequency response node that strongly attenuates a sinusoid, the total power of narrow-band noise centered at the same frequency will not be attenuated as much. It will survive with a dip in its spectrum. Room modes or frequency response peaks are usually not as severe and narrow as frequency response nodes.

An example speaker and room magnitude frequency response (CC-BY-SA 3.0 by Gregory Maxwell)

However, I wonder if your application affords the degradation of audio quality due to the modulation. Also, the room frequency response is not affected so if the test audio is noise, you won't get any sort of improvement.


The impulse response, $h(t)$, of a system I have this differential equation
  $$
    y''(t) + y'(t) + y(t) = x(t) - 3x'(t)
    $$

I found the system response of the system:

$$
        H(f) = \frac{1 - 6\pi fj}{(2\pi fj)^2 +2\pi fj+1}
    $$
which is equivalent to

$$
        H(s) = \frac{1 - 3s}{s^2 +s+1}
$$

My question is: How to find the impulse response? If I factorize the denominator, it's not gonna be pretty. e.g., $s + \frac{1 +\sqrt{-3}}{2}$

Also, what does the following mean?

Find the output for $t \ge 0$ with initial conditions $y(0) = 1$, $y'(0) = -1$ when the input is $2cos(6t)u(t)$. Identify the zero state, the zero input, the transient, and the steady state responses.

firstly do partial fraction for H(s) and then use transform tables of Laplace Transform. Using inverse laplace and tables you must be able to find impulse response. I think you can use j^2 instead of -1 and find complex roots for denominator. Then, using transform pairs you will find a complex solution.


How to perform a time domain shift in the frequency domain without zero padding Apologies in advance if my question doesn't make sense. I'm not the most fluent person in signal processing.

I have an array of time domain values that I would like to apply a time delay to by computing the fft then apply the Fourier transform
$f(t - \tau) = F(\omega)\exp(-i \omega\tau)$ and then compute the inverse fft to get my delayed signal.

My problem is that the delay is very large (three times further than the range of my original signal) so in order to get the delayed signal I have been zero padding my original signal by a significant amount which greatly increase the computation time.

Is there a better way to compute the delay using Fourier transform without the need to zero pad?

You must zero-pad, whether implementing the delay in the time domain or the frequency domain. (Consider this: by delaying, you are making the signal longer.)

Implementing the delay with the FFT implements a circular shift. If you don't pad and you use the FFT, the data will simply wrap around on itself. (Imagine that if you didn't zero pad and used the FFT with $\tau = NT$, the duration of the record length. Then you would simply get back what you started with, with no delay, because of circular wrapping.) If you want only an integer number of sample delays, do it in the time domain.

If you want a fractional sample delay, you can use the FFT as you describe. It makes no difference mathematically if you implement say a 20.7 sample delay first in time by 20 samples and then in frequency by 0.7 samples, or whether you do it all in frequency by 20.7 samples. Remember—you padded first. Computationally, as you say, this can increase computation time. But you are only increasing the length by a factor of four so computation time with an FFT should increase by a factor of only two. Is this too much? Alternatively, you can do the "bulk" delay in time and the fractional delay in frequency. In the frequency part, you have a couple of choices. First, before doing the bulk time delay, you can do the fractional delay on the original unpadded data and accept the fractional-sample delay wrap-around error at the beginning of the sequence. Second, you can pad your data by only one sample and do the fractional shift with an FFT. If your FFT software accepts only powers-of-two record lengths and your original record length is a power of two, this will likely slow your computation by more than a factor of two because a fast algorithm is not employed. However, most modern FFT packages provide fast algorithms for sequence lengths that can be factored into products of small prime numbers, so adding one data point might not significantly increase computation time.

Some have suggested alternatively using an interpolater for a fractional shift. The FFT is the ideal interpolater—it exactly interpolates bandlimited data when zero padding without adding a delay exponential or when your delay exponential is added. However, cubic splines are outstanding interpolators and should be in your toolbox.


Low samples numbers machine vision training for object detection I'm trying to set-up a machine vision application in which I have to detect the position of a printing nozzle (in the red sqaure in the image) in a field of view in order to prcisely position a ganty stage.

I currently struggle to detect the wished features. I have tried to implement a Haar cascade using opencv but the detection is not performing.

To train my classifier I used the full image with the location of the nozzles blacked out as negative and the cutout of the red rectangle as positives.
I have few positives and negatives (15 of each) that are basically the same image with different illuminations/expositions and placements of the printer head in the field of view.

Since this approach did not work, what kind of image recognition should I use in that case? Would a training of a custom set in tensorFlow give better results? I've seen implementations of DNN using existing vast libraries of tagged images but I'm not sure how if this is the best method for this application.

Thanks in advance

If the camera has a shallow depth of focus and the background is mostly blurry and the foreground is only moving at one focal plane, why not use a simpler method, such as Template Matching?

Since the printing head is changing its angle of view with respect to the lens as it moves in the $x,y$ dimensions, you might have to use more than one templates but overall the process is straightforward: Apply the cross correlation in a sliding window across the image and stop where the cross correlation peaks.

This example is in MATLAB but easily transferable to other platforms.

This sliding window might seem slow but is unavoidable. Neural networks do it implicitly at the lowest level of inputs. The good news here is that the printing head image does not move in the $z$ dimension, in which case the scale would change too. This would mean that you would have to repeat this process for the scales of interest.

You might be able to derive "hints" of where the printing head is roughly in the image, by tracking its silhouette. If you threshold the image for example, you might want to look for a distinct pattern of large rectangular areas in black, that roughly match the size and shape of the arm. This would help in making the search area where the normalised cross correlation is applied, smaller. But overall, I would not expect a modern Desktop or even SBC to struggle with an HD image and a patch even as big as 128x128 pixels.

Hope this helps

EDIT:

In terms of a pesudo-code instruction:

Take a representative image of the nozzle (don't crop it too much around the nozzle but make sure the bit you want to recognise is clearly visible).
Given an image from the monitoring camera, run a normalised cross correlation across the image. How you are going to do this depends on the platform you are using but MATLAB, Octave, Python all have means of running a normalised cross correlation, fast, in this way
The output of this will be a Real two dimensional array that will have a single peak exacty at the point where the template had the  best match.
Threshold that match and take the centroid of the region (if it spans more than one pixel) to find out where the nozzle is in the camera image.


Is interpolation (interp1) better than FIR filtering when rational integers are close to 1? Question
I've been attempting to resample a GPS signal in MATLAB. I've built a few FIR filters using fvatool and from handmade transfer functions (punched out with a HP-35s). Most are kaiser-windowed LPFs, but some are least-squared LPFs. All of the filters I've built have nominal responses that should adequately filter the aliasing out of the interpolations. I've also built a FFT-Resampling tool that resamples well and very quickly. However, I've hit a strange anomoly, at least as I see it. A simple 'nearest' neighbor interpolation using interp1 seems to much more accurately resample my signal than any FFT-Resampling I have done, or any upsample, filter, decimate tool I've used, to include dsp.FIRRateConverter, resample,upfirdn, and even straight upsample, filter,decimate. Why is this? Can anyone explain why this occurs? My theory is that when the rational integers of your sampling ratio are very close to one, straight interpolation will provide more accurate results because you reduce occilations and/or other noise distortion.

Background
My signal is a standard collected GPS signal, pulled from the air, sampled at 25e6 sps and saved in a binary file using signed 16-bit integers. When pulled into MATLAB (has to be done in chunks due to size) the signal is a complex row vector.

I'm upsampling to 26.25e6 sps, therefore my rational integers p and q are 21 and 20, respectively. To use interp1 I simply build two time vectors:

 t0 = (0:1:numSamples - 1)*1/sampleRate;
 t1 = (0:1:numResamples - 1)*1/resampleRate;

where

 sampleRate= 25e6;
 resampleRate = 26.25e6;
 samplingRatio = resampleRate/sampleRate;
 numSamples = length(s0) % length of my original signal, sampled at 25e6 samples
 numResamples = length(s0)*samplingRatio % length of my sampled signal; should equal length(s1)
 s1 = interp1(t0,s0,t1,'nearest','extrap');

When using a FIR filter, I have used many methods. But the simpliest is to call resample like this (resample is a least squares FIR filter with a kaiser window called through firls that gets pushed through upfirdn):

 s1 = resample(s0,p,q);

where

 [p,q] = rat(resampleRate/sampleRate,1e-12);

You are looking at the wrong metric of "correctness".  Nearest neighbor is introducing significant discontinuities that are showing up as massive quantities of noise in the result.

The problem is that you should be comparing to the result you would have gotten if you had sampled at 26.25MHz in the first place.

Let's try it (sample a 12Hz sine wave at 26 Hz for 2 seconds):

numResamples = 52;
resampleRate = 26;
t1 = (0:1:numResamples - 1) * 1/resampleRate;
s1correct = sin(t1 * 24 * pi);
plot(t1,s1correct);
plot(t1,abs(fft(s1correct)));

Ok.  Now let's see what happened if we had sampled a 12Hz signal at 25Hz, and then interpolate up to 26Hz using matlab's linear interpolation.  (Actually, I'm using Octave because I'm on a machine without matlab today)

numSamples = 50;
sampleRate = 25;
t0 = (0:1:numSamples - 1) * 1/sampleRate;
s0 = sin(t0 * 24 * pi);
s1linear = interp1(t0, s0, t1, 'linear', 'extrap');
plot(t1,s1linear);
plot(t1,abs(fft(s1linear)));

Pretty good.  Now let's look at nearest neighbor interpolation.

s1nearest = interp1(t0, s0, t1, 'nearest', 'extrap');
plot(t1,s1nearest);
plot(t1,abs(fft(s1nearest)));

It's completely different, and filled with noise at all frequencies!!!  (The "beat pattern" in the nearest neighbor result is the beat pattern from sampling at 25Hz instead of 26Hz.)  The noise is because nearest neighbor is compressing the 12Hz signal to 12.48Hz, and then inserting an extra sample every 26th sample to get us back to 24 peaks.  That extra sample is the discontinuity that is creating all that extra noise.

I think you'll find that cubic interpolation is even better than linear.


Property of the trace and expectation I'm reading the paper Model-Driven Deep Learning for Joint MIMO Channel Estimation and Signal Detection by Hengtao He, Chao-Kai Wen, Shi Jin, and Geoffrey Ye Li on Orthogonal Approximate Message Passing (OAMP) detectors.
There is an step that it is getting me crazy, in Equation $(42)$ in Appendix B. It is the following:

At first, I thought that it was a consequence of a Hypothesis, but then, to derive the next expression in Equation $(44)$, it is used again this formula.

Is there some kind of property which I might be missing?
Regards.

I am not really sure which part you didn't understand. I'm assuming it's the step that goes from the rightmost expression at the top to the one at the bottom.
Assuming $\mathbf{q} = \mathbf{x} - \hat{\mathbf{x}}$ is i.i.d.
Let me give it a try...
$
\begin{aligned}
\textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) ]
& = \textrm{tr} \{ \textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) ] \} \\
& =  \textrm{E}[ \textrm{tr} \{(\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) \} ]  \\
& =  \textrm{E}[ \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \} ]  \\
& =   \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} ~\textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}}) (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}} ] \}   \\
& =   \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} ~\textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}}) (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}} ] \}   \\
& =   \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} ~\sigma_{q}\mathbf{I}] \} \\
& =   \sigma_{q} \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H}] \}.
\end{aligned}$
where I am assuming $\sigma_{q} = \textrm{E}[  \|(\mathbf{x} - \hat{\mathbf{x}}) \|^2 ] = \frac{\textrm{E}[  \|\mathbf{q} \|^2 ]}{N_t} $.
As you can see, I am making a lot of assumptions here. I don't know if they hold since I couldn't access the paper. But I think it can shed some light on your problem (if your question is strictly about the trace operations).


How to solve this Laplace integral for an averaging function in an iterative numerical (finite difference) model? I am trying to implement a mathematical model for vibrational damping described in this article.
They provide an equation for damping force ($F$) as a function of:

a spring constant ($k$)
a damping coefficient ($η$)
displacement of the point being damped ($x$)
velocity of that point ($\dot{x}$)
mean displacement and velocity for the point ($x_m$ and $\dot{x}_m$)

I understand the basic principle of the force equation. However, my math is not very advanced, and I don't understand how they calculate the mean displacement and velocity ($x_m$ and $\dot{x}_m$), given by the equation for $z_m$.
The relevant excerpt is as follows:

The equation for $z_m$ (where $z = x$, $z= \dot{x}$) is what I don't understand. I presume $n$ here is a time constant that allows a relative weighting for older or newer values. Is that correct?
What does the expression $t^{n+1}$ represent? Do they mean time raised to exponent of $(n+1)$? Or something else like time at sample $(n+1)$?
What about the integral? Do they mean a Laplace $s$? What do I do with this in a stepwise function? They explain that this ought to be solved using a stepwise iterative approach which is perfect because I am using stepwise finite difference modeling for the system I want to add this to.
In a finite difference model where $Δt$ is the time step per time sample, $n$ is a given coefficient for mean weighting, and where you have the displacement $x$/$x_m$ and velocity $\dot{x}$/$\dot{x}_m$ at the prior time step, in plain English or simple demonstrative code or line by line math, what would be an iterative solution they are suggesting for $z_m$ (ie. the current $x_m$ or $\dot{x}_m$)?
Thanks for any help. It is very much appreciated.

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


LTI system and initial conditions I suspect this is rather obvious once explained; I must be misunderstanding something.

From reading various answers here (such as this one), I understand that an LTI system must have zero output for zero input, which I think is equivalent to zero initial conditions or initial rest.  To quote @MattL from the answer I referenced,

  "This explains why a system with non-zero initial conditions can
  neither be linear nor time-invariant. A linear system must have a zero
  output for zero input"

In the textbook I am working through, one of the problems states "Find the response for the system described by the following difference equation with the initial conditions given", and the difference equation is as follows:

$$y[n]-0.7y[n-1]=u[n], y[-1]=-3$$

Here the input starts at $n=0$ yet $y[-1]=-3$.  The textbook is asking the question in the context of LTI systems.

At first I thought there was a mistake in the question, but I have seen similar setups in answers here on stackexchange as well.

How is the above scenario compatible with zero initial conditions?

The system is described by a linear difference equation with constant coefficients and as such, it is described in the same way as a linear time-invariant system. It is just the non-zero initial condition that makes the system non-linear, at least according to the common definition of linearity in system theory (homogeneity and additivity).

Such a system with non-zero initial conditions is also referred to as incrementally linear. An incrementally linear system responds linearly to changes in the input. For this reason it is very common, and it also makes a lot of sense, to discuss it in the same chapter as general discrete-time LTI systems.

As a final note, in practice we're mainly interested in stable systems. For such systems, the influence of non-zero initial conditions becomes negligible after some time, and, consequently, the system will practically behave as an LTI system after the transient caused by the non-zero initial conditions has died out.


Harmonic errors in acoustic instruments If a note has a fundamental frequency of 100Hz, then its second harmonic is located at 200Hz and its tenth harmonic is located at 1000Hz. But for acoustic instruments, the harmonics will not occur at exact multiples of the fundamental. My question is, how much error is typical and what is the relationship between the degree of error and the harmonic number. So for example, is it common for the error in the tenth harmonic to be as much as, say, 50Hz. 1050Hz in other words. That's a 5% error. But would a 5% error be just as likely for the second harmonic?

A few things:

Not all instruments have partials neatly defined by the $h_n = nf$. Kettle drums would be an obvious example, but bells are another. In this case, you will be way off looking in the "usual places" +/- 5%.
Many instruments have only even or odd harmonics. Others are designed to omit some harmonics. The piano, for example, is designed to have almost no 7th harmonic. This may not matter to you, but it's worth pointing out in case it does.
It is common for higher harmonics to be "out of tune", that is, not in the "usual places". An engineer who once worked as a piano tuner once admonished me for using harmonics to tune my upright bass (a technique that works well for lighter-stringed instruments), because, he said, the partials are more out of tune with heavier stringed instruments. He said it's especially true told of higher harmonics, explaining that the math that gives you the simple formula for partials assumes that strings have no mass. 5% strikes me as a little much (almost a semi-tone!), but I don't think it's out of the question.
Many instruments have vibrato. This throws all your harmonics out of wack:


Effect of window function on audio output If I have a low pass filter, does the choice of window function I use have any effect on the quality of the audio input? What would be the effective difference between using a rectangular window and a blackman window of the same size?

A cursory search of google will give you frequency response plots of various windows, for example: http://www.labbookpages.co.uk/audio/firWindowing.html


IIR Lattice Filter in C programming I would like to implement Lattice IIR filter in c, i can't defined the value of gN, i know that x(n) = fN(n), but how to initialize gN ?. you can find the algorithm here : n

assume all of this (except k[0] and the output y[]) gets defined.  states g[] are initialized to zero.
int N, NUM_SAMPLES;
float g[N+1], v[N+1], k[N+1];
float x[NUM_SAMPLES], y[NUM_SAMPLES];

then the sample processing would look like this:
int n = 0;
while (n < NUM_SAMPLES)
    {
    float f = x[n];

    float acc = 0.0;
    int m = N;
    while (m > 0)
        {
        f -= k[m]*g[m-1];
        g[m] = k[m]*f + g[m-1];
        acc += v[m]*g[m];
        m--;
        }
    g[0] = f;
    acc += v[0]*g[0];

    y[n] = acc;
    n++;
    }

that is, i think, the essential C code.


A question about the sampling theorem Let's say I have a signal, containing frequencies up to $30$ kHZ.

I am sampling this signal at $40$ kHz with infinite bit depth (just for theory).

Will I be able to recover exactly all the frequencies up to $20$ kHZ , or is that the signal must  only contain frequencies below $20$ kHZ to get the exact magnitudes of the frequencies up to $20$ kHZ ? In other words, will frequencies above $20$ kHZ introduce uncertainty about those below $20$ kHZ ?

If you sample your signal with $f_s=40\,\text{kHz}$ then the band between 20kHz and 30kHz will be folded back to the band between 10kHz and 20kHz. So aliasing will affect the upper half of your spectrum after sampling, i.e. the band between 10 and 20kHz. The lower band from DC to 10kHz will not be corrupted by aliasing.


what are the Pros and Cons of ESPRIT and MUSIC algorythms? On wikipedia there are only 2 lines about ESPRIT: Estimation of signal parameters via rotational invariant techniques (ESPRIT)... MUSIC doesnt have an entry in the music disambiguation page, other than something for neural networks.

I have read that they are more advanced than FFT. What are their advantages and disadvantages?

The prime thing such algorithms aim to do is to make use of more information that you may have about the signal. In this case, the extra information is that you know the number of signals (sinusoids) present in your measurements.

One pro for both is, therefore, when your measurements match the assumption, you get a more accurate representation of the spectrum: it has the right number of lines. The FFT does not guarantee this.

A con is that, if your assumption is wrong (there are three lines instead of two), then these algorithms perform worse.

These algorithms do not give better resolution than the FFT. Only more data can give more resolution.

There is a good exposition here comparing the two to each other.


designing FIR narrowband pass filter using windowing technique I m designing FIR narrowband pass filter using windowing technique. is bandpass filter is same as narrowband pass filter? I m not finding any good article on this.

At 44.1 kHz sampling frequency, an FIR (and possibly also linear phase) digital filter with a pass band of 20 Hz would require very large number of taps. If this is not an issue for you then the design is straight forward by the windowing method. You can use Kaiser window design principle to find the minimum length of a Kaiser window that meets your design criteria (namely the maximum pass and stop band ripples, and the transition band width).

The following OCTAVE / MATLAB commands design a very narrow band bandpass FIR filter,specified by its parameters according to your only given cirteria.
The filter length is not unique, you can play with it to find the minimum
that meets your other criteria (if there are any).

Fs = 44100;        % Sampling Frequency
fc1 = 490;         % cutoff frequency 1
fc2 = 510;         % cutoff frequency 2

L = 20000;         % length of the FIR filter
h = fir1(L , 2*[fc1,fc2]/Fs ,'bandpass');    % bandpass FIR filter.

figure,freqz(h,1,L);
title('Frequency response of the FIR filter');

and the result is: (I've zoomed into the pass band frequency region to observe it clearly...)


carrier frequency difference Papers use wireless device-dependent radio-metrics as fingerprints. A radio-metric is a component of radio signal, amplitude, frequency and bandwidth. Each device creates a unique set of radio-metrics in its emitted signal due to hardware variability during the production of the antennas, power amplifiers, ADC and DAC circuits.As a result, radio-metrics cannot be altered post-production, and thus provide a reliable means for distinguishing wireless devices.

In fact, the transmitted RF signal from a wireless device experiences hardware impairments, channel characteristics, and noise at the receiver. Distortions caused by channel-specific and noise-related effects are likely to have a more random structure. Distortions in a feature that are caused by transmitter hardware impairments should manifest themselves consistently across multiple frames from the same transmitter.

device-dependent fingerprints are Modulation based features like carrier Frequency offset, phase shift offset, I/Q origin offset in the constellation plane.

How can I simulate carrier frequency offset in MATLAB?
In the following reference 1, the carrier frequency difference (CFD) is used as a feature for distinguishing wireless devices and it is assumed as Gaussian distribution. Why can we assume that this feature is Gaussian distributed?
Why do noise and the channel not affect this feature?

On Identifying Primary User Emulation Attacks in Cognitive Radio Systems Using Nonparametric Bayesian Classification

//Edit

CFO of different devices in the feature space

//reference

Device fingerprinting to enhance wireless security using nonparametric Bayesian method

1 As defined, the device-dependent fingerprint offset should be constant for a given device so if you have nominal_carrier_freq = 1000; then  you could say

unique_dev_offset = 1;
unique_dev_carrier_freq = nominal_carrier freq + unique_dev_offset

2 Using the Gaussian distribution is common for modeling any deviation from nominal. Think of the context. Our goal is to have devices that are identical, but in real life they have some slight differences. Take a resistor for example, with a nominal value R. If we have a set of resistors we will see some are a little different than R but the majority of them are centered around R. So our model single modal distribution around R. With decreasing probability of being far away from R. This could be a lot of functions (triangle for instance meets those requirements) but the Gaussian is used because it has nice mathematical properties and usually fits real world data very well

3 I think you have to say the channel is constant for this to work. now imagine we have two signals. both are supposed to be the same but they have some unique characteristics. in the images below, I have two signals, one at nominal frequency (0 carrier freq offset), and one with a nonzero carrier frequency offset. We can see that even after adding noise the signals are still distinguishable from each other. Noise won't affect the characteristic frequency offset of the signals. It may modify the appearance of the signals, but not the underlying characteristics

hope that helps


Computing SSD/SSE using pattern with omitted pixel The goal is to compute a weighted SSD over a small neighborhood of pixels (Np pattern in the following figure)

I want to know why the bottom-right pixel is omitted:

It is said that to enable SSE-optimized processing but I didn't really understand it!

This is the link for the PDF: http://vision.in.tum.de/_media/spezial/bib/engel2016dso.pdf

and it is the figure - 4

By excluding the bottom right pixel, the authors only need to process 8 pixels which can be stored in eight SSE registers (xmm0 through xmm7).
https://en.wikipedia.org/wiki/Streaming_SIMD_Extensions#Registers


Discrete Filter $y[n] = \frac{1}{3} x[n] + \frac{1}{3} x[n-1] + \frac{1}{3} x[n-2]$ Consider the filter which equation can be represented by $y[n] = \frac{1}{3}x[n] + \frac{1}{3}x[n-1] + \frac{1}{3}x[n-2]$, in $x[n]$ and $y[n]$ are sequence of input and output of the system respectively.

Mark the correct option about the sequence ${h[n]}$ filter response with respect the unit pulse ($\delta_0[n] = 1$) in $n=0$.

A) h[1] = 2/3

B) h[1] + h[5] = 1/3

C) h[0] = 0

D) h[2] = 3

E) h[2] + h[3] = 1

My attempt:

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

$1+z^{-1} + z^{-2}$ its a Finite Filter Response right?

$X(z) = 1$

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

I stuck here because and really cant realize to envolve the discrete domain with this expression on Z domain.

Any hints?

You can see by definition that your system is LTI System.
Moreover it is casual as it is only depends on input from the past.

The quaestion asks for response on the Unit Pulse signal on time 0.
In LTI Casual System it is the coefficient which responds to to input at time $ n $ -> Its value is $ \frac{1}{3} $.

I'm not sure where $ h \left[ 1 \right] + h \left[ 5 \right] $ is from though.
As clearly this system can be represented by 3 Coefficients which you found in Z Domain (Though they are also given from the finite differences equation).


Best digital noise reduction in optical spectra measurement I am using an optical spectrometer to measure the reflectance of some surfaces; the output is essentially a set of photodiodes' signals: for each photodiode I have a discrete time signal every $T$. The exposure time and the frame rate of the sensor can be changed: I tuned them so to have the highest number of measurements per second (i.e. highest framerate possible) and the corresponding longest integration window possible (there is no saturation problem).

Given as a degree of freedom the number of measurements that can be acquired, what would be the best digital method to reduce the noise of the measurement and to maximize the SNR?

The simplest solution seems to be averaging, which enhances the SNR of $\sqrt N$ under hypothesis of white uncorrelated noise, but is that the best that can be done in case of unknown noise spectral power distribution? In particular, I was wondering if there was a technique that exploits in an optimal way the information gathered on the noise during the acquisition of $M$ repeated samples.

For the additive white noise channel, it can be shown that matched filtering maximizes SNR; and for constant illumination, the "transmitting pulse filter" is a rectangle, and the match to that is a rectangular weighting, i.e. an average, too.

The idea is simply that you "my signal is a vector; multiplication / convolution then becomes a dot product. Cauchy-Schwarz Inequality gives me an upper bound for the norm of signal·receive window / filter, and I reach that when what I do at the receiving end is a linear multiple of what I do at the transmitting end".

For unknown spectral distributions, you'd try to find a filter on the receiving end of things that is as orthogonal as possible to the spectrum of the noise. You'd need to estimate that first, though! Devise some calibration method; but, really, when you're observing the same spot for long enough, maybe just look at the spectrum and ignore the constant value?

Assuming you have some non-white noise, you can "abuse" techniques from the world of spread-spectrum communications design. Imagine the following:

You can modulate the intensity of your light source to two values, a and b, over time. Imagine sending a pseudo-random sequence consisting of times when you send the higher intensity, and of the lower intensity:

a      ----    --------    ----

(a+b)/2

b          ----        ----

At the receiver, the first thing you do is simply subtract the average of the received signal intensity; that way, what you receive becomes:

c      ----    --------    ----

0

-c         ----        ----

Now, that's a very simple spread spectrum signal (in fact, DSSS), and it's almost certain that its spectral properties are pretty different from these of noise.

Now, you just correlate with the high-low sequence that you've sent, considering "high" a "1", and "low" a "-1". When you shift that "+1 -1 +1 +1 -1 +1" sequence correctly, you'll see that by point-wise multiplication between "c -c +c +c -c +c" (which you received), you get a very high amplitude: 1·c + (-1)·(-c) + 1·c + 1·c + (-1)·(c) + 1·c = 6c . In reality, you use lengths much higher than 6!

You still get the same gain as when averaging white noise over a length of $N=6$, but for low-frequency-heavy noise distributions, you'd get noise suppression, too. (you buy that by getting less noise suppression at higher frequencies! No such thing as a free lunch; but you might want to do that, because you might have other light sources in your room – for example, 100 Hz modulating room light, which you can eliminate efficiently).

A bit of general commentary:

That averaging gives you the gain you state depends on the assumption that noise is uncorrelated with signal! In general, that's not the case for you: The energy of noise is very much dependent on current forward current in a semiconductor junction; that's kind of a bummer. What the spread-spectrum approach does is in fact a bit of noise shaping, and it might well suppress a bit of the harmonics that nonlinearity introduces, but I'm not expert enough in photonics to tell you how severe these nonlinearities are.

Also, I don't know the reflectivity of what you're measuring. Just be warned that white LEDs aren't really white: they consist of multiple light sources that generate different wavelengths; same goes for your photo sensors: They aren't uniformly sensitive to all wavelengths. If you haven't already, calibrate your sensor with a color proof sheet!


How do I find Gray Level Variance of the Image? Tell me how do I find out gray level variance of the image? Is it any different from the normal variance of the image.

I would think that it is indeed the "normal" variance of the image.  You look at the distribution of the pixel values (i. e. gray levels), and compute its variance.


FFT of a periodic function Suppose you have N samples of a sine wave of frequency $f_0$ sampled at $f_s$.

$$
x[n] = A\sin[2\pi(f_0 /f_s)n]
$$

How do you get a delta function in the frequency domain (i.e. all of the energy is in the center one bin of the FFT)?
Does the number of samples (N) need to be a prime number of cycles of $f_0$? If so, why?
Does $f_s$ need to be an integer multiple of $f_0$? If so, why?

In reverse order,

Yes the sampling frequency needs to be a multiple of $f_0$.
The number of cycles of the sinusoid in $N$ samples needs to be an
integer, not necessarily a prime number.
If you do 3. and 2., all the energy will be concentrated in two bins
(corresponding to $100$ Hz and $-100$ Hz.

Be careful, though. If you have a sinusoid $x(t)=\sin(2\pi (100) t)$ at 100 Hz
and you sample it at 1000 Hz, you will have $10$ samples per cycle of the
waveform. With $x[0] = 0$ and $x[9] \neq 0$, you can do an "FFT"
of length $10$ if you like, and all the energy will be in $X[1]$ and $X[9]$.
I put FFT in quotation marks because at such short lengths, a plain
vanilla brute-force DFT calculation might be almost as fast. Note also
that $x[10] = 0$ is not part of the DFT calculation.
Alternatively, you can use $N = 10n$ samples and do an FFT of length $1000$  (say) on samples $x[0]$ through $x[999]$.  In this case, the
energy will all be in $X[10]$ and $X[1014]$.
Regardless of what some other answers will likely tell you,

  DO NOT zero-pad the length $1000$ FFT to length $1024$ to get an
  "efficient" implementation.

If you do follow the advice on zero-padding, your length-$1024$ FFT
will NOT have the energy concentrated in two bins as you desire.


OFDM synchronization using training symbols In OFDM receivers, special OFDM training symbols for which the data content is known, can be used to perform the synchronization process. A correlation operation is performed between the received signal and the training symbols. From the correlation peaks, both symbol timing and frequency offset can be determined. I have two equestions.

Practically, how these OFDM symbols are generated at the receiver?
If these reference symbols are sufficient for both symbol tiiming and frequency offset determination, why differential detection is often employed?
The information bits are modulatedon the subcarriers by differential QPSK. The reference phase and amplitude of each subcarrier have to be known in order to estimate the bits correctly. For this purpose differential detection is used.

Claim

In OFDM receivers, special OFDM training symbols for which the data content is known, can be used to perform the synchronization process.

Be a bit careful about that statement.
There's basically four kinds of synchronization that I can think of right now:

Frequency synchronization
Timing synchronization
Phase synchronization
Clock synchronization

... as well as higher-level things such as frame sync, decoder sync.
Frequency sync
Can be done in a number of ways. A very OFDM-typical way of doing so are autocorrelation methods.
With these methods, you don't actually use a known reference symbol, but send OFDM symbols that have a specific temporal structure.
For example, the well-known Schmidl&Cox synchronization is based on the idea that if you have two preamble symbols with certain symmetry, you can do both a fine frequency estimate as well as a timing estimate.
Rough frequency estimates are often done by exploiting the fact that OFDM systems usually don't use all subcarriers, but leave guard bands at the edges. So, for this, too, no known symbol is needed.
Timing sync
... to expand on the above: It's very easy to build a timing estimator based on a fixed lag correlation, when you have a symbol that repeats itself. You just compare $N$ samples with $N$ samples a fixed delay $\Delta T$ later; when that correlation function sees a peak, you might be detecting your preamble symbol. Side effect is that the phase of that correlation (divided by $\Delta T$) serves as fine frequency estimate.
DAB (and DAB+), for example, use a non-Schmidl&Cox autocorrelation based frequency and timing estimator. WiFi uses S&C, for example.
Clock Sync
That's a tough one with OFDM, as it's usually pretty bad if the receiver sampling clock is different from the transmitter's clock. It's solved by limited symbol length, or by coupling the carrier frequency to the sampling clock (so that you can correct the sampling clock by correcting the frequency), or by estimates based on the phase sync.
Phase sync
This is really the only thing that is normally done through known data. But it's commonly not done by known symbols, but only by known pilots within distributed "ressource blocks", i.e. single subcarriers within payload symbols, as that has a lot of advantages by observing the channel as its being used.

Practically, how these OFDM symbols are generated at the receiver?

not at all. They are received, just like every other OFDM symbol. That's why it's called a receiver.
The known pilot subcarrier symbols are typically written in some standard document.

If these reference symbols are sufficient for both symbol timing and frequency offset determination,

They aren't. Frequency synchronization is usually done with autocorrelation methods, see above.

why differential detection is often employed?

Are they? Unbacked claim! I know that DAB does use OFDM with differential PSK internally, but that is kind of an odd duck: It's from the early 90s, so that OFDM was feasible on small systems, but proper channel estimates weren't. Generally, using OFDM says "need to use a wideband channel efficiently", and that usually means that you really want that non-differential modulation BER goodness.

The information bits are modulatedon the subcarriers by differential QPSK.

You're referring to some specific system, not to something that applies to OFDM in general. So, this is, in general, a false statement.

The reference phase and amplitude of each subcarrier have to be known in order to estimate the bits correctly. For this purpose differential detection is used.

This is self-contradictory. Either you need to know the phase of each subcarrier, or you do differential modulation. Never both.


Why are FFT numbers so high? My input numbers are close to 10, here is a sample of some:

9.259992,9.253408,9.322236,9.375503,9.406026,9.507771,9.509567,9.403034,9.397048,9.36054,9.265977,9.262985,9.216302,9.228272,9.338395,9.430565

I am using this PHP library
thinkingmik/FastFourierTransformation

My peak numbers in response are about 500. A control which was done using another tool is 1-1.5

FFT is well over my head, so I am trying to read this library to get a bit of an idea, but can't figure this out. The interesting thing is the client says the chart seems to be correct, just numerical values are wrong.

As pointed out by Dan Boschen in a comment, these number aren't really surprising. Usually the DFT is defined as

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}\tag{1}$$

where $x[n]$ is the input sequence and $N$ is the DFT length. So bin zero, i.e., $X[0]$, is just the average value of your input sequence multiplied by the DFT length. The maximum value of $|X[k]|$ can be as large as the maximum value of $|x[n]|$ multiplied by the DFT length.


what sampling rate and bit depth should i use? I'm doing a simple mini project on FPGA board, electronic lock which open when a short track (the "password") is play next to the lock's mic.

i want to save samples of the chosen song to an array x[n]
and and then take samples from the mic and save them to a FIFO list y[n].
then compare the the conv result z[n]=x[n]*y[n] and when it is big enough the lock is open.

the question is how many samples i need for 3-4 sec of rock song for this purpose?

The typical audio sampling rate is 44 kHz, so you would need between 130 000 and 175 000 samples.
In practice though, it depends on the sampling rate you use on your FPGA. You can lower, but you will need to low-pass your signal before sampling it to avoid aliasing.

Note that if you put a threshold on z[n]=x[n]*y[n], it will not be invariant to volume differences.


Z transform convergence From the definition I know:

$$|\sum_{n=-\infty}^{+\infty}x[n] \cdot z^{-n}| < \infty\space\space\space(1)$$

Is there any sequence of x[n] which can not be written as

$$x[n] = y[n] \cdot u[an-b]\space,\space\space\space a,b \in \mathbb{Z}$$
(where the dot means a multiplication and $u[k]$ means the discrete step function. In other words, x[n] is an infinity sequence)

but do satisfy the definition in (1)?

Nop. For an infinite sequence, the convergence of Z transform depends on the value of $z$ also. So as far as $n\rightarrow\pm\infty$, equation (1) will not be satisfied for all values of $z$.

Consider a simple case, $x(n) = u(n)$
$$X(z) = \sum_{n=0}^{\infty}z^{-n}$$
$$X(z) = 1 + z^{-1}+z^{-2}+z^{-3}+ \cdots$$
for $z=0.5$,
$$X(z) = 1 + 2+4+8+ \cdots = \infty$$
for this sequence to converge, $|z|>1$.


How to make sense of audio file's $\textrm{dBFS}$? I'm trying to build a song visualizer using the FFT method. When I analyze the audio, I get the $\textrm{dBFS}$ value for each frequency range.

From my understanding,

$\textrm{dBFS}$ is in the range $[-\infty, 0]$, where $0$ is the most loud.
The formula for percentage loudness is $10^{\rm dBFS/20}$.

However, the $\textrm{dBFS}$ values I'm getting are consistently below $-20$, i.e., $10\%$ loudness, even at loud parts of the song. Why is this so? Am I calculating the loudness wrongly?

Also, how can I calculate the overall loudness at a point in time? Is it the sum of all the values, or the maximum?

However, the dbFS values I'm getting are consistently below -20, i.e., 10% loudness, even at loud parts of the song. Why is this so? Am I calculating the loudness wrongly?

Why should that be a problem? Human perception is roughly logarithmic. Sound that isn't "unbearably loud" shouldn't be much more than 10% of your "acoustic dynamic range".

By the way, with $\frac{}{20}$, this formula gives you amplitude, not power values, so this isn't truly loudness, just amplitude. But since that, in perception, is but a constant factor, these terms can be used interchangably.

  Also, how can I calculate the overall loudness at a point in time? Is it the sum of all the values, or the maximum?

It's the sum, at least approximately: human perception has a different sensitivity for different frequencies, and also, frequencies can mask each other, so the energy in the signal is only proportional to the perceived loudness if we ignore these facts (see: psychoacoustic models).


Extended Kalman Filter (EKF) Linearization of Non Linear Functions I have more general question about the extended Kalman filter usage. What is not clear to me is why the EKF uses non-linear functions $f$ and $h$ for state prediction and estimate, while in other places the Jacobian of these functions is used.

Why the following is never used?

First calculate the linearized state and measurement models at previous estimate point using Jacobian. Use the linearized state transition and measurements matrix everywhere instead of non-linear in this specific iteration.

We can't use the initial linearized-model of non-linear system again and again in each iteration because linear-model approximates the nonlinear model only around some small neighbourhood of the starting point. In other words, linear model is good(i-e approximates the nonlinear model) only around some small-interval around the point. As, we need to estimate the current state again and again, if we don't linearlize the non-linear model around the current estimate in each iteration, we would we getting worse and worse (in terms of linear model) and consequently we will no longer correctly estimate the state.


Why are some LDPC codes more immune to noise than others? When one analyzes an LDPC code, it is usual to plot the BER versus the $E_b / N_0$. My questions are:

Why are some LDPC codes better than others in terms of noise immunity? After all, they are just defined by matrices $\mathbf H$.
What properties or parameters of this parity-check matrix have influence on the BER for a given SNR?

As you mention, the LDPC code is completely determined by its generator matrix H. Hence, the properties of H define the (theoretic) performance of the code.

Since the LDPC code is a linear code, the ultimate measure for the (theoretic) code performance is its minimum distance between two codewords, that is:

$$
d = \min_{c_1,c_2\in\mathcal{C}}\|c_1-c_2\|_0
$$

where $\mathcal{C}$ is the set of existing codewords (i.e. bit patterns that fulfill the parity check) and $\|\ldots\|_0$ denotes the hamming distance. This is equal to the expression (since it's a linear code)

$$
d = \min_{c\in\mathcal{C}}\|c\|_0
$$

where $\|\|_0$ is the Hamming weight.

I wrote the word "theoretically" in parenthesis, since the finally obtained performance does also depend on the decoding algorithm. Due to the low-density of the parity check nodes, the LDPC is suitable for iterative decoding, which can yield big gains.

Another improtant measure is the length of the shortest cycle of the LDPC Tanner Graph. The shorter the cycle, the worse is the performance of the iterative decoding.


Discretized CWT vs. DWT vs. FWT vs. MRA By reading some books about wavelets, it seems that Discretized CWT is not the same thing as DWT.

Can we give a classification of Discretized CWT vs DWT vs FWT etc. ?

A few remarks :

All I read about DWT or MRA (multi-resolution analysis) involves a dyadic factor (2)  for $a$ in $\frac{1}{\sqrt{a}} \Psi(\frac{t-b}{a})$, ie $a=2^m$. Is it always true ?
Remark : This is good in order to have a non-redondant transform (see B Z's answer here Scalogram (and related nomenclatures) for DWT?), but the graphical representation of such a transform isn't very satisfying for visual analysis of a signal.
The fact that there is $ a=2^m$ involved implies that the frequency bands studied [in the case where the mother-wavelet has a narrow-band fourier transform, e.g. Morlet wavelet http://ieeexplore.ieee.org/ieee_pilot/articles/06/ttg2009061375/assets/img/article_1/fig_3/large.gif] are things like [20hz,40hz], [40hz,80hz], [80hz, 160hz], etc.
Thus we cannot hope a good graphical representation with precise frequency resolution. It seems that this kind of transform is more suited for "approximation + details" than for visual frequency analysis (scalogram, etc.)
On the other hand, a discretization of Continuous wavelet transform allows to take $a=1.0004^m$ for example, and thus have a more precise frequency resolution, for one who wants to do a graphical representation (scalogram).

For the CWT the values in your definition can be any $a\in \mathbb R^+$ and $b\in \mathbb R$ if applied to functions. The CWT can also be applied to discrete signals (discretized CWT) with $b$ restricted to appropriate discrete values (and any $a$, but typically $2^{j/v}$ where $v$ is the "voices/octave" of the transformation). Now when we restrict $a$ to some discrete values $2^j$ we get the DWT, which is an orthogonal transformation.

(Small edit to correct the scale : $2^{j/v}$ instead of $a^{j/v}$)


What is the meaning of BPSK(n) I am doing some research into GNSS (Global Navigation Satellite Systems) and have come across various signals listed as having the in-phase or quadrature channels modulated with BPSK(2) or BPSK(10), etc.  However, I am struggling to find a definition of BPSK that includes this additional parameter.

What does the (n) signify in BPSK(n)?

I understand the basic concept of BPSK - Binary Phase Shift Keying and have implemented this in very basic tests by multiplying the LO by 1 or -1 for a 1 or 0 chip respectively.

It has been suggested by a colleague that the n is simply the chip rate and that BPSK(1) is a chip rate of 1MHz, BPSK(2) 2MHz, BPSK(10) 10MHz, etc. as mentioned in the comments by Jim Clay.

While I do not have proof, it seems pretty clear, given the "Characteristics of Compass signals reported as of May 2008 compared to GPS-L1CA" table at the Wikipedia page that you offered as an example, that the "(n)" notation refers to the chip rate.  You will notice that there is one BPSK(1) signal that has a chip rate of 1.023 Mchips/s.  There are four BPSK(2) signals listed, and all of them have a chip rate of 2.046 Mchips/s, which is twice 1.023 Mchips/s.  There are also three BPSK(10) signals that have a chip rate of 10.23 Mchips/s which is, of course, 10 times 1.023 Mchips/s.

It appears that the "(n)" notation means that the signal has a chip rate of $n*f_o$, and in the case of the Compass signals $f_o$ is 1.023 Mchips/s.


How to represent a given equation more clear, profressional and short form? I have a image $U_{m \times n}:\Omega \to \mathbb R^2$, the output $P$ can be define as
 $$P=\mu J_{m \times n} - U$$
where $\mu = \max \{ u_{ij} : 1 \leq i \leq m, 1 \leq j \leq n\}$, $J_{m \times n}$ be the $m \times n$ matrix whose $i, j$th component is $1$: that is, the all-ones matrix. (This notation isn't quite standard, but it's as close to standard as I know. $J$ is often the all-ones matrix)

However, it is so many sentences for expression the above equation. Do we have more short and standard way to represent it? As my found, the $J_{m \times n}$ can be expressed by indicator function such as

$$P=\max (U)\times 1_{\Omega}-U$$

where $1_{\Omega}$ is Indicator function

Does it equivalent with original meaning? If not, please give me a standard and common expression in image processing. Thanks

I think the notation $J_{m,n}$ is pretty standard for an all-ones matrix. So the only possible clarification would be to dispense with $\mu$ and directly specify what it actually is:

$$P=||U||_{\max}J_{m,n}-U$$

where $||.||_{\max}$ denotes the max norm of a matrix.


How do I find a system's impulse response from its state-space repersentation using the state transition matrix? Suppose we have a linear represented in the standard state space notation:

$$ \dot{x}(t)=Ax(t)+Bu(t)$$
$$y(t) = Cx(t) + Du(t)$$

In order to get its impulse response, it is possible to take its Laplace transform to get

$$sX=AX+BU$$
$$Y=CX+DU$$

and then solve for the transfer function which is

$$\frac{Y}{U}=C(sI-A)^{-1}B+D$$

Similarly, for a discrete system, the $\mathcal{Z}$-transform of
$$ x[n+1]=Ax[n]+Bu[n]$$
$$y[n] = Cx[n] + Du[n]$$

is

$$\frac{Y}{U}=C(zI-A)^{-1}B+D$$

This process seems a bit long and I remember that there is a way to find the impulse response using the state transition matrix which is the solution for $x$ of the first equations of each pair. Does anyone know how to do this?

You can approach the problem using the state transition matrix by solving the standard non-homogeneous ODE in the first equation. The solution to $\dot{x}(t)=A x(t) + B u(t)$ is

$$x(t)=x_0 e^{At}+\int_{0}^te^{A(t-t')}Bu(t')dt'$$

where $x_0=x(0)$. The quantity $e^{At}$ is called the state transition matrix (also the solution to the homogeneous ODE), which I'll refer to as $\Xi(t)$ (I don't recall the standard notation for this). Taking $x_0=0$, the equation for $y(t)$ becomes

$$y(t)=C\int_0^t\Xi(t-t')Bu(t')dt'+Du(t)$$

The above equation gives you the output as the input convolved with the system impulse response and indeed, you can take the Laplace transform of the above equation to verify. Noting that the Laplace transform of $\Xi(t)=e^{At}$ is $(sI-A)^{-1}$ and convolutions in the time domain become products in the s-domain, we get

$$Y=C(sI-A)^{-1}BU+DU$$

which gives you the same transfer function as in your question.

Regarding your comment on the fully Laplace transform approach being long, I wouldn't necessarily say it is so. However, state transition matrix approach might be simpler to implement, because several operations involving it can be computed with simple matrix multiplications and nothing more.


camera calibration change orientation of axis I want to solve for the extrinsics by using direct linear transformation on corresponding 3D LIDAR points and 2D camera points. I already have the intrinsics.

Problem is, points behind the camera gets re-projected as well (see picture below).

So I constrain to only points "in front of the camera", i.e z > 0. The problem is, on different trials where different sets of points are used, the produced extrinsic matrix produces differing axes. Sometimes, constraining z > 0 gives the right results (centre part of image), whereas other times I need z < 0. So the question is, how do I constrain the Z axes of the camera to be sticking out of the camera?

def with_intrinsic(points2d, points3d, intrinsic):
    cam1_K_inverse = np.linalg.inv(intrinsic)

    #direct linear transformation calibration, assumes no intrinsic matrix
    assert points2d.shape[0] >= 3
    assert points3d.shape[0] == points2d.shape[0]

    A = []

    points2d_homo = []
    for u,v in points2d:
        points2d_homo.append([u, v, 1])

    points2d_homo = np.array(points2d_homo).T #columns to be data points

    points2d_inv = np.dot(cam1_K_inverse, points2d_homo).T
    assert points2d_inv.shape == (points2d.shape[0], 3)
    assert points2d_inv[0, 2] == 1

    for idx in range(points2d.shape[0]):
        x3d, y3d, z3d = points3d[idx]
        u, v, _ = points2d_inv[idx]

        A.append([x3d, y3d, z3d, 1, 0, 0, 0, 0, -u * x3d, -u * y3d, -u * z3d, -u])
        A.append([0, 0, 0, 0, x3d, y3d, z3d, 1, -v * x3d, -v * y3d, -v * z3d, -v])

    A = np.array(A)

    U, D, VT = np.linalg.svd(A)
    M = VT.T[:, -1].reshape((3, 4))

    error = get_reprojection_error(points2d, points3d, intrinsic, M)
    logging.debug("error with_intrinsic: %s", error)

    return M

The question is a bit unclear (and am sorry to say, the code even more so). Assuming your 3d-2d point correspondences are correct, all you'd need to do after solving the DLT is to rearrange the orientation of the axes while leaving their directions and the parity of the frame unchanged (i.e. ensure that the determinant of the rotation matrix is +1).


Problem with Convolution in MIMO System For a flat fading channel (SISO system) we have:
$y=h \star s+n$

But for a mimo system, we have:
$y=Hx+n$ with $H \in \mathbb{C}^{N_r \times N_t}$

Why is there not the product of convolution for a mimo system?

I think it should be:

$$ \boldsymbol{y}_{i} = {H}_{i} \boldsymbol{x} + \boldsymbol{n}_{i} $$

For the $ i $ -th antenna and $ {H}_{i} $ being the convolution matrix of this specific channel.
You could write this in a Matrix Form where $ H $ becomes a tensor.
In case $ \forall i, j, \; {H}_{i} = {H}_{j} $ then it can be made simpler.

You could also work on the samples:

$$ \boldsymbol{y} \left[ m \right] = H \left[ m \right] \boldsymbol{x} \left[ m \right] + n $$

Where in this case $ {y}_{i} \left[ n \right] $ is the $ i $ -th antenna at time sample $ m $ and the $ i $ -th row of $ H \left[ m \right] $ is the response of the channel $ i $ at time $ m $.
This model allows time variant channels.


How we can find ROC(Region of convergence) given a signal in z domain? I read that ROC is a region,which is a set of values where z transform is defined ,that is it converges

Lets say i have a discrete time time signal $x[n]=n^2 u(n)$ and i want to find its ROC(Region of convergence),how can i do that?

You need to determine the values of the complex variable $z$ for which the series

$$\sum_{n=1}^{\infty}n^2z^{-n}\tag{1}$$

converges. So you need to know a few things about infinite series.

For this case a simple test is the ratio test. You take the ratio of two successive elements of the series and compute the limit:

$$L=\lim_{n\to\infty}\left|\frac{(n+1)^2z^{-(n+1)}}{n^2z^{-n}}\right|=\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\left|z^{-1}\right|=\frac{1}{|z|}\tag{2}$$

The series converges absolutely for $L<1$, i.e., for $1/|z|<1$, or $|z|>1$. So the ROC is the region $|z|>1$.

Note that the ROC would not change if we used any other power of $n$ in $(1)$.

EDIT:

A step for step explanation of Eq. $(2)$:

$$\begin{align}L&=\lim_{n\to\infty}\left|\frac{(n+1)^2z^{-(n+1)}}{n^2z^{-n}}\right|\\&=\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\left|\frac{z^{-(n+1)}}{z^{-n}}\right|\\&=\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\left|z^{-1}\right|\\&=\frac{1}{|z|}\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\\&=\frac{1}{|z|}\lim_{n\to\infty}\left(\frac{n^2+2n+1}{n^2}\right)\\&=\frac{1}{|z|}\lim_{n\to\infty}\left(1+\frac{2}{n}+\frac{1}{n^2}\right)\\&=\frac{1}{|z|}\left(1+\underbrace{\lim_{n\to\infty}\frac{2}{n}}_{0}+\underbrace{\lim_{n\to\infty}\frac{1}{n^2}}_{0}\right)\\&=\frac{1}{|z|}\end{align}$$


A better explaination for 'there's no flipping of the impulse response' I showed the question about flipping the impulse response to some math teachers and electrical engineers at an university. They read Dilip's answer and thought that it was hard to understand what Dilip meant and it was simply not true...

Could someone confirm that it really isn't true what Dilip Sarwate wrote. And explain it in a way that is easier to read.

See: Flipping the impulse response in convolution

  There is no "flipping" of the impulse response by a linear (time-invariant) system. The output of a linear time-invariant system is the sum of scaled and time-delayed versions of the impulse response, not the "flipped" impulse response.

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


Which algorithm could be used to compare two waveforms from an oscilloscope in XY mode? I would like to receive some hint as to which algorithm could be used to compare two waveforms from an oscilloscope in XY mode, as shown in the image:

Note that the image shows a percentage difference (displayed in red) between the two curves (one yellow, one blue).
It is this percentage value that I would like to obtain.
I understand that the curves are produced point to point (x, y) and then a sequence of position data (x, y) is generated:

I believe that the algorithm would have to compare two data arrays.
Does anyone have any idea how to do this?
Thank you.
Edit: Following the formula of the video (Normalized Correlation), I believe that the result of comparison is interesting, my next step is to try to do it with a microcontroller to find out if it will have resources.
Note: After calculating the correlation for each axis, I multiplied the result of the axes and inverted the result, to show the difference. If this is not correct, please comment on any suggestions.

Ref.:

Normalized cross-correlation is also the comparison of two time
series, but using a different scoring result. Instead of simple
cross-correlation, it can compare metrics with different value ranges.

Source: https://anomaly.io/understand-auto-cross-correlation-normalized-shift/index.html
Edit (03/14/2021):
In addition, the measurement is performed with a dedicated microcontroller and synchronized with the frequency of the sine wave, that is, the number of readings is fixed.
The basic circuit:

The name “Octopus” maybe doesn’t say anything alone but if you google
it along the words “curve tracer” you will obtain thousands of
result.In few words an “Octopus” curve tracer is a small circuit that
used in conjunction with a scope allows to display  the voltage across
a component under test on horizontal (X) axis versus the current
through that component on the vertical (Y) axis.A scope set to X-Y
mode is required (most of them have this feature).

This   circuit will produce a “signature” waveform on the oscilloscope
to aid in the testing and analysis of shorts, opens, and leakage in
just about any electronic component including resistors, capacitors,
inductors, diodes,transistors and digital ICs too.Each component has a
characteristic waveform (called “Lissajous” pattern), some examples:

Source: https://www.jammarcade.net/simple-component-tester-a-k-a-octopus-curve-tracer/

You can't get the exact value because you don't know how log the signals spent at each stage.
Consider the case where you have two square with period 1, and shifted by $\delta$ compared to each other.
The XY plot will always be a rectangle with the width corresponding to the amplitude of the X signal and height corresponding to the amplitude of the Y signal. In other words, you have always the same input for your function to compute the correlation, while the correct output varies.
You still can get an idea of how similar two signals are, but it is not the signal cross correlation


Orthogonality of FSK modulation for non-coherent detection For an FSK modulation with coherent detection, in order to ensure orthogonality, the frequency separation between symbols is $ \Delta f=f_2 - f_1 =1/2T$. However, in case of non-coherent demodulation, which must be this value $1/2T$ or $1/T$?

For continuous-phase FSK, frequencies $f_1$ and $f_2$ and phases
can be chosen so that a frequency spacing of $\frac{1}{2T}$ gives orthogonality.
Note that it is not sufficient to have just $\Delta f = \frac{1}{2T}$; both
$f_1$ and $f_2$ have to be integer multiples of $\frac{1}{2T}$ or $\frac{1}{4T}$ etc.
Look for Sunde's FSK and minimum-shift keying (MSK) regarded as coherent FSK.
If the phases are not carefully specified (as they indeed are for Sunde's
FSK and MSK), orthogonality requires $\Delta f$ to be an integer multiple
of $\frac{1}{T}$ and so of course the minimum spacing is  $\frac{1}{T}$.

The theory of noncoherent detection assumes that the phase of the received
signal is arbitrary rather than carefully controlled, and it is guaranteed
to give the specified error-rate vs SNR performance regardless of signal phase
as long as the frequency-spacing is $\frac{1}{T}$. The noncoherent receiver
will work  and it will provide the specified error-rate vs SNR performance
if its input is a continuous-phase FSK signal with minimum frequency spacing $\frac{1}{2T}$ etc. However, it is worth keeping in mind that if the received
signal is such that coherent detection is possible, then use of noncoherent
detection in lieu of coherent detection is a suboptimal choice as far as BER
vs SNR is concerned, though it might be preferable for other reasons such as
ease of implementation, cost, weight, power consumption, etc of the receiver.


Audio Steganography: Inaudible Audio Watermarking for Source Identification What is a relatively inaudible audio watermarking technique to overlay a broadband audio signal on top of playing music without analysing the music for low-bitrate source identification? Ideally should be resistant to multi-path, reverb and resonance.

I'm not trying to do any kind of DRM, so I don't care about resampling, time compression or frequency shifting - I just want to figure out in real-time which source is playing without adding too much audible noise to the signal.

At the moment I'm broadcasting a low bitrate ultrasonic signal on top of the music, modulated with GMSK around a carrier frequency of 19500Hz, but the high frequencies attenuate fast and the high frequencies are still audible.

I'm considering broadcasting a super wide-band PN-sequence that sounds like pink noise. I only need to sustain a transfer rate of about 3-4 bytes (24-32 bits) per second. But if I stretch out the symbol length, I'm worried about environmental noise.

I could add a delay if I need to pre-analyse the output buffer before, but I hope that I can introduce a very soft full-spectrum "hiss" or hum that can carry my source identification signal.

I don't need to track the source during periods of silence, so there will always be some music to mask my signal.

Update: I've found some prior audio steganography work, but I don't know what the trade-offs are for real-time modulation (esp. on a low-power embedded device) and if this would at all work without look-ahead on the output buffer.

Petrus. First of all, thanks for attaching my github library as an update into your question for other readers. I have just seen your question.

According to my own experiences, echo hiding methods are fragile towards additional echoes, so are towards reverb. A simple spread-spectrum model using a PN Sequence, which introduces some noise, embeds data into high-frequency bands, so hidden data is usually audible and not robust against low-pass filters.

Transformations such as cosine transforms and wavelet transforms are often used in audio watermarking. Applying a spread-spectrum method without a PN Sequence onto Approximation Coefficients of the third level wavelet decomposition (A3) would be very robust against many common attacks except high-pass filter, and would be almost completely imperceptible. On the other hand, hiding data onto Detail Coefficients of different levels, gives different balance between robustness & transparency. These decompositions are done very easily in Matlab if you have Wavelet Toolbox. However, I have not uploaded my related works into github yet.

We have not worked with real-time audios, but I know that all of these techniques can be applied for them as well. I am not sure whether you already solved your problem, but if you have any further questions you can point out in a comment. I am not an expert in digital signals unlike many of you on this site, so different opinions are also very welcomed.

Best Regards.


Designing a Notch filter for Digital Audio to be able to work with all audio sampling frequencies I am currently designing an FIR notch filter using Matlab Filter Design & Analysis tool. I created a notch filter at desired notch frequency (lets say f0) with sampling frequency of 44.1kHz. But I found that the same notch filter is not able to filter-out/attenuate the desired frequency in other audio signals with different sampling rates.

Is it possible to design a single notch filter that works for all audio sampling frequencies.

FYI, I have tried with designing a notch filter with highest sampling frequency, but it is not working for audio with sampling rate other than what it is designed for?

Can somebody help me/provide me a lead in designing such a notch filter?

What you want is impossible. The notch frequency is always relative to the sampling frequency, so you design a notch filter with a given ratio $f_{notch}/f_s$. What you need to do is design different notch filters for different sampling frequencies. The same is of course true for cut-off frequencies of frequency selective filters, center frequencies, etc.

By the way, considering computational complexity, IIR filters are usually better for implementing notch filters than FIR filters.


Camera Calibration - why do you have to move the calibration board? I understand the mathematical derivation of Zhang's Method for camera calibration but what I don't understand is why we need to move and rotate the calibration board for several frames in order to get accurate calibration results? One obvious reason is that in order to detect radial distortion you want the control points to be uniformly distributed along the field of view of the camera. However, couldn't we have also just ensured that the calibration board fills the FOV of the camera and just keep it static? Why do we have to move it around?

Thanks!

This is due to the optimization problem being rather high-dimensional (around 11 parameters). With only a single observation of the calibration board, there would be multiple possible combinations of parameters explaining the observed feature point locations (unless a very constrained camera-model is used). Only a sufficient number of sufficiently independent observations will properly constrain the optimization to yield a well-defined answer. Otherwise the routine will converge to some (local) minimum.

Think of this as solving a set of equations with more unknowns than equations. You will be able to find a solution, however, there would be an infinite number of solutions (one or more free parameters).

With respect to a board viewed from the front (frontoparallel orientation), filling the entire image. Now -- is this a calibration board close to the camera viewed with a wide-angle lens (short focal length), or is it a calibration board far away, viewed with a zoom lens (long focal length)? This shows that we need images with the board tilted (foreshortening).

With some constraints, you can calibrate a camera with a single image of two calibration patterns forming e.g. a 90 degree angle, or two perfectly parallel boards at differents depths.

I have written an article on the usual optimization process, which might help you understand how the process could end up not converging to the global solution: https://calib.io/blogs/knowledge-base/camera-calibration

This has little to do with the fixed grid of the imaging sensor as suggested in other answers.


What do different Cepstral or MFCC coefficients represent intuitively? I understand the explanation for separating slow and rapid changing log spectral components but i need to understand:

Why lower coefficients have higher and mostly positive magnitudes?
Why Higher coefficients have  negative values more
often?

I understand the mathematical computation i.e. taking log of spectrum to make the components with different change rates, linearly seperable and then taking IFFT of logspectrum to separate these components.
However, I analyzed the MFCC stats from a speech dataset and can't see why slowly changing components (lower coefficients) have more magnitude while higher ones have more negative values and what do negative valued coefficients indicate.

When you take the log of some measurement value, you are simply changing the nature of the scale of what you are measuring.  Where the dividing line between positive and negative values on the new scale is depends only where the "1" was defined on the original scale.  There is no qualitative distinction between positive and negative per se, so what you are really asking is why are slowly changing components (lower coefficients) larger than other ones.
If a sinusoidal is enlarged proportionally, then the amplitude is inversely proportional to the frequency.  If your signal indicates this (regardless of context) it means similarity is being preserved across levels.
So this doesn't answer the title question, but hopefully helps with the enclosed ones.


Maximum Throughput of QAM system Consider a QAM system designed to transmit over a bandwidth of $3 \ \mathrm{kHz}$. The channel's power constraint imposes a maximum $\mathrm{SNR}$ of $30 \ \mathrm{dB}$. The system can tolerate a probability of error of $10^{−6}$. I want to determine the maximum throughput of the system in bits per second.

Let's call the throughput $T$. If $M$ is the constellation size and $"$ is the bandwidth, then I know that $T=MW$. I've calculated the value of $M$ using theformula $$M = \log_2 \left( 1 - \frac{1.5\mathrm{SNR}}{\ln(p_e)}\right)$$ Then I approximated $M$ to the nearest power of $2$, which is $8$.
Therefore, $T = 24000 \ \mathrm{\frac{bits}{sec}}$
but it is a wrong answer. Can anyone tell me what I have done wrong?

When $\log_2(M)$ is even, the relationship betweeen BER and QAM if an AWGN channel is given as below (Assuming Gray coding so that a symbol error in most cases is one bit error):

\begin{align}
k&=\log_2(M)\\
y &= 10^{\frac{SNR}{10}}\\
P_e&=\frac{4}{k}\frac{\sqrt{M}-1}{2\sqrt{M}}\textrm{erfc}\left(\frac{\sqrt{3k\frac{y}{M-1}}}{\sqrt{2}}\right)
\end{align}

Using the above formula to solve for $y$ given your SNR and constellation $M$ should achieve the result you are looking for.

This is an excellent explanation for the derivation of the formula, specific to QAM.

And here is a nice summary by the same author for BER vs SNR of many different modulation schemes


Higher-order moment of output of LTI system Assume a very simple LTI system. Assume $x$ is white Gaussian i.i.d. with variance $\sigma^2$.

The output variance is straightforward to obtain. For example, for a continuous-time system:

$$\mbox{Var}(y) = \mathbb{E}(y^2) = \sigma^2 \int_0^{\infty}|H(s)|^2 \,\mathrm d s$$

However, I am interested in higher order moments, for example, $\mathbb{E}(y^4)$. How can this be obtained (is there a closed form solution in the first place)?

PS: Note that $\mathbb{E}(y^4)\neq\sigma^4 \int_0^{\infty}|H(s)|^4 \, \mathrm d s$.

For the case of input process $\{X(t)\}$ being white Gaussian noise with two-sided power spectral density $\frac{N_0}{2}$, the output process $\{Y(t)\}$ is a strictly stationary zero-mean Gaussian process in which all the random variables have the same variance $\frac{N_0}{2}\int_{-\infty}^\infty |H(f)|^2 \,\mathrm df$ almost as you say.  But the key point here is that the $Y$'s are zero-mean Gaussian random variables with known variance, and it is known that if $Y \sim N(0,\sigma_Y^2)$ then
$$E[Y^n] = \begin{cases}\sigma_Y^n
\cdot (n-1)
\times  (n-3)
\times (n-5)\times  \cdots \times 3 \times 1, & n ~ \text{even},\\
0, & n ~ \text{odd}.\end{cases}$$

Now, if the input process is just white noise (not necessarily Gaussian white noise), then the output process is not necessarily Gaussian, and while the mean and variance are as specified above, we cannot say anything about the distribution of the $Y(t)$'s and cannot infer anything about the higher moments of the $Y(t)$'s the way we could for Gaussian output processes.

If the input process is a zero-mean WSS Gaussian process (which is also a strictly stationary process) but not necessarily a _white noise process, than $\{Y(t)\}$ is also a zero-mean WSS Gaussian process (and so strictly stationary too) with variance
$$\sigma_Y^2 = \int_{-\infty}^\infty S_X(f)|H)f)|^2 \,\mathrm df$$ where $S_X(f)$ is the power spectral density of the input process and so once again we are in business and can use $$E[Y^n] = \begin{cases}\sigma_Y^n
\cdot (n-1)
\times  (n-3)
\times (n-5)\times  \cdots \times 3 \times 1, & n ~ \text{even},\\
0, & n ~ \text{odd}.\end{cases}$$


What is the difference between frame error rate (FER) and symbol error rate (SER)? What is the difference between frame error rate (FER) and symbol error rate (SER)?
I know the difference between frame error rate (FER) and bit error rate (BER) from the website below

https://www.quora.com/What-is-the-difference-between-frame-error-rate-FER-and-bit-error-rate-BER-in-LDPC-codes

In this website,it said

Bit is the unit that data transmit in the physical line. Frame is consist of many bits.
A unit frame has more information than a unit bit.

However,the symbol consist of many bits too. e.g.symbol "0" = 000,symbol "7"=111,they all consist of three bits,so what is the difference between frame error rate (FER) and symbol error rate (SER)?

Let's say that you transmit a signal using 4-PAM, with symbols $3,1,-1,-3$. The bits-to-symbols mapping is

$\begin{align*} 3 &\rightarrow 00 \\ 1 &\rightarrow 01 \\ -1 &\rightarrow 11 \\ -3 &\rightarrow 10.\end{align*}$

The symbol error rate (SER) is the probability that the receiver estimates the wrong symbol; for example, that it estimates $3$ when the transmitted symbol is $1$. Note that this is different than the bit error rate, because some symbol errors cause one bit error, while others cause two errors (for example, estimating $3$ when $-1$ was transmitted).

This is relevant for the physical layer in a communications system. However, in the data link layer (DLL), the "protocol data unit" is a frame; that is, the DLL can only deal with entire frames. A single bit error in an entire frame counts as a frame error; usually, the frame is discarded and a retransmission is requested. The error is typically detected by a CRC. The frame error rate is the ratio of errored frames to transmitted frames.

As an example, consider communication using BPSK, and a channel with probability of bit error $P_b = 1\times10^{-3}$. The DLL uses 64-bit frames. The probability that a frame has zero errors is $$(1-P_b)^{64} \approx 0.938,$$ which means that the probability of receiving an errored frame is $$\text{FER} \approx 1-0.938 = 0.062,$$ or one errored frame out of every 16 transmitted frames.

You can reduce the probability of handling an errored frame to the DLL using forward error correction in the physical layer. If you use a Hamming (7,4) FEC code, which can correct one single error in a 7-bit word, then the probability of bit error out of the decoder is reduced to $$P_d \approx \binom{7}{2} (1-P_b)^5 P_b^2 = 2.1 \times 10^{-5}.$$

Using the same 64-bit frames, the probability that a frame has zero errors is $$(1-P_d)^{64} \approx 0.9987,$$ which means that the FER is $$\text{FER} \approx 1-0.9987 = 0.00067.$$ In other words, one frame out of every 1500 will have at least one bit error.


How to calculate signal which is not changed by a filter? Suppose that there is a FIR filter F and a signal S.  The filtered signal is the convolution of F and S, F * S.

The problem: how to calculate a signal S' such that F * S' = S' (the filtered version is the same as itself, up to a fixed time shift), and S' minimizes norm(S, S') (where norm is some kind of distance metric, for example sum of squared differences).

Loosely: how to calculate the most similar signal to a given signal which is not affected by a filter?

EDIT: The original problem may be over constrained because there may be very few (or no) signals for which F*S = S exactly.  An alternate formulation is to find S' which minimizes the combined error alpha*norm( F*S' - S') + beta*norm( S' - S ), for some weights alpha and beta.

EDIT: The specific filter I have in this case is F=[ 0.00097656, -0.00976562,  0.06347656, -0.15625   ,  0.18554688,   0.83203125,  0.18554688, -0.15625   ,  0.06347656, -0.00976562,    0.00097656] which is a fairly modest lowpass filter.  I was interested in the general case as well though.

If a filter with transfer function $H(f)$ passes a signal $x(t)$ with Fourier transform $X(f)$ unchanged, then $H(f)X(f) = X(f)$ for all $f$, and
so $H(f)$ must have the property that
$$H(f) = 1 ~ \text{for all frequencies}~ f~ \text{for which}~ X(f) \neq 0.$$
Thus, an arbitrary filter will not enjoy this property for very many
signals, perhaps for no signal at all.  Indeed, if $H(f_0) = H(-f_0) = 1$ for
some $f_0$, then the filter will pass the sinusoid
$A\cos(2\pi f_0t + \theta), -\infty < t < \infty$ unchanged, but it will
not pass a sinusoidal pulse such as
$$x(t) = \begin{cases} A\cos(2\pi f_0t + \theta),
& -T < t < T,\\ 0, &\text{otherwise,}\end{cases}$$
without distorting it. Of course, it is possible that $H(f) \neq 1$ for
all $f$ and so even passing a monochromatic sinusoid is not possible
(cf. @hotpaw's comment).
In short, the subspace referred to by @AtulIngle might well be
the trivial subspace of zero dimension containing only the
zero signal.
It is also worth noting that if $H(f)$ has constant nonzero value
over an interval, then the filter is physically unrealizable for
essentially the same reason that an ideal low-pass filter is
unrealizable.

So, the OP's question is likely to have no answer at all because the
signals $S^\prime$ that he seeks to find may not exist at all for
the given filter. Perhaps something else was intended to be asked?