Symbol timing synchronization using a high sampling rate Symbol timing synchronization usually relies on a digital PLL procedure, an iterative process that relies on detecting a timing error signal that is used to control an interpolation filter that in its turn tracks this error until lock is achieved.
 Suppose this feedback algorithm is totally ignored and brute force oversampling is used so that the maximum amplitude sample at the output of the matched filter would be the one to process further. so no tracking, no lock. Ignoring any small quantization error, are there any other issues using this method?

The primary issue with this approach is that it will appear to work fine in your design if you do not consider noise, but in actual application will fail under low or even moderate SNR conditions unless you start doing what an appropriate tracking loop response will do by design and that is average out the sample to sample noise to provide your best estimate of the correct timing location; that which you would not be able to properly estimate on any one given sample. So yes as long as you do the proper filtering/averaging and timing estimation and correct your timing estimate based on current and past updates you will be able to do your timing recovery- but this is what the timing recovery loop is doing and very efficiently.  For any such design that you pursue for timing (and carrier recovery, and data demodulation) be sure to test and benchmark your approach using corner case SNR conditions AND worst case symbol clock offset between the transmitter clock, signal as received, and receiver clock!

To add to this since not detailed in my other links is the Gardner Timing Error Detector (TED) implementation suitable for QPSK and QAM implementations, where the symbol $y_n$ shown is the complex I,Q output. I have not implemented this for OQPSK specifically, but do not see why simply offsetting the I and Q samples half a symbol would not work. The TED produces an average level that is proportional to timing error and will have significant pattern noise on any one given output so is assumed to be averaged as part of the overall loop design:

Note the the Gardner TED has a relatively high capture range for carrier offsets (see Location of Matched Filter for more details on that), so assuming the estimation of the carrier is within this range, the timing recovery can be done ahead of carrier recovery as shown in this QAM receiver implementation (two samples per symbol for timing recovery, which once acquired is decimated to one sample per symbol for carrier recovery. Also depending on the waveform pulse shaping used and timing error determination the inputs to the timing recovery may be better before or after matched filtering (hence the dashed lines in the block diagram showing that either may be used).

Below is a graphic that compares an approach of changing the frequency of the sampling clock to the ADC in the receiver vs interpolation and resampling. The latter is more prevalent in modern software radio solutions.

(It is unfortunate that the two implementations are in reverse direction flow and I apologize for that; it is the figures I had readily available)

The following graphic demonstrates the concept of interpolation and resampling :

A polyphase filter approach as publicized by fred harris would do this interpolation and retiming efficiently, so that you never actually need to interpolate to a higher rate than 2 samples per symbol. Note that you also only need three filter implementations in the poly-phase filter bank: a present and the two adjacent. After changing from present to adjacent, the filter coefficients can be updated from a ROM table. (The reason for keeping the three is so that all filters needed are properly loaded with the waveform history prior to switching two its output).

For further details on polyphase retiming and implementation of polyphase filters see these links:

Emulating a Variable Delay

How to implement Polyphase filter?


Zero mean noise in images According to the literature, zero mean noises are those whose expected value is zero. But what does this mean in the context of images. Let us consider grayscale for now. The intensity values ranging from 0 to 255. The noise will also have values between 0 and 255. How can then the mean be zero?

To get a nice visualisation lets consider a colour image with zero mean noise & standard deviation 20! Which means the values of each RGB component can randomly vary from +20 to -20. So pixels will be added with this noise.

Here is the code for generating noise to an Image using opencv c++

Mat mSource_Bgr,mNoise_Bgr;
mSource_Bgr= imread(FileName_S,1);

// We need to work with signed images (as noise can be
// negative as well as positive). We use 16 bit signed
// images as otherwise we would lose precision.
Mat mGaussNoise_Bgr(mSource_Bgr.size(), CV_16SC3);

double StdDev=10;
randn(mGaussNoise_Bgr, Scalar::all(0), Scalar::all(StdDev));
Mat mTemp;
mSource_Bgr.convertTo(mTemp,CV_16SC3);
addWeighted(mTemp, 1.0, mGaussNoise_Bgr, 1.0, 0.0, mTemp);
mTemp.convertTo(mNoise_Bgr,mSource_Bgr.type());

mGaussNoise_Bgr.convertTo(mGaussNoise_Bgr,CV_8UC3);
imshow("Noise Image",mGaussNoise_Bgr);
imshow("Original Image",mSource_Bgr);
imshow("Original Image + Noise",mNoise_Bgr);

Here is the Original Image

Noise Image with Mean : 0 & Std Deviation : 20

Original Image + Noise

Same Image with Noise Mean:0 Std Deviation:10


Best Guess Rectangle I am using openCV to find a particular object, and am trying to find a way to make it locate the object consistently. The object's location will be determined by taking 30 frames and trying to find the best guess based on where it is found in each of the 30 frames. Both the camera and object are (mostly) stationary (the object can move slightly, but optimally it will be stationary during the 30 frames).

This image demonstrates what I'm looking for. I know that the green box may not truly be what would be calculated, I just threw it in there to demonstrate. The red boxes are the guesses on each frame, and the green one is what I'm trying to calculate. Are there any known formulas/algorithms for this, or any references on where to start?

Note

I initially tried to just calculate the mode along the X-axis (ignoring width) (with each estimate rounded so that there would be some similarity), and picking a random rectangle whose x value was in that range for X but that was not giving me a very good estimate.

I must admit I haven't worked on anything like this before, but looking at your data it seems to me that computing the median of e.g. the lower left corner point of all rectangles (i.e. taking x and y dimensions into account) should get rid of outliers. After having removed outliers from your data, you could simply try to average the remaining rectangles. Obviously, such an approach requires some heuristic tweaking, but it may be worth a try.

I just tried Matt's approach (just taking the 1D medians of the top/left/bottom/right X and Y values) and the results seem pretty good. See plot below of 20 realizations of various random rectangles.

scilab CODE ONLY BELOW

function r = getRect(t,l,b,r)
    large =  [1 1 0 0]*rand(1,1,'normal')*abs(t-b) + [0 0 1 1]*rand(1,1,'normal')*abs(l-r);
    small =  [1 1 0 0]*rand(1,1,'normal')*5 + [0 0 1 1]*rand(1,1,'normal')*5;

    noise = small;
    if (rand(1,1,'uniform') < 0.2) then
        noise = large;
    end
    r = [t l b r] + noise;
endfunction

for p = 1:20;

N = 20;

x_points = [2 2 4 4 2];
y_points = [1 3 3 1 1];

subplot(5,4,p)
RR= [];
for n=1:N
    R = getRect(10,10,100,50);
    RR = [RR; R];
    plot(R(x_points),R(y_points),'r');
end

RRm = median(RR,'r');

xset("thickness",4)
plot(RRm(x_points),RRm(y_points),'g');
xset("thickness",1)
end;


Find the Equivalent Noise Bandwidth I have already asked (in another thread) how to quantify (specify) the window function I have created. One figure I need it seems is the Equivalent Noise Bandwidth (ENBW) Or Noise Equivalent Bandwidth. I now need some practical help. How would I go about determining the ENBW of my window function? I guess I am asking, what signal I should apply and what should I measure to generate some numbers from which I can then calculate the ENBW? Having read around the subject, I have not found a very straightforward explanation of what to do. Apologies for being really dense, mathematics is tough for me because of my form of dyslexia. Thanks.
This first image shows my window (black trace) for which I would like to calculate the ENBW. This is shown versus Blackman Harris (green trace) window for reference. This test emulates a perfect (synthesised) 1kHz sine wave sampled at 48kHz with 32-bit precision and 1Hz bin-width.

This second image shows my window again (black trace) versus a Blackman Harris again (green trace) only with an attempt to emulate a realistic situation. The synthetic sine wave is now 1000.5Hz (so 47.976 wavelengths fit in one sampled block) and a noise floor has been added. This shows how my window allows the noise floor (and second tone - see below) to be viewed in detail while it is obscured by the Blackman Harris window. Also, a second sine wave of 1050.5Hz and -160dB has been added to demonstrate the point of the improved window - to better separate components that are close in frequency. Please note that levels are all over the place because of a lack of calibration rather than level errors in the window itself.

SUMMARY
The equivalent noise bandwidth (ENBW) for a window function is the bandwidth in bins of a brickwall filter that would result in the same noise power from a white noise source as the DFT "filter" (when viewing, appropriately, each bin of the DFT as a bandpass filter). The ENBW for the rectangular window (no further windowing) is 1 bin as demonstrated in the first plot below. The ENBW for any window can be determined from the following equation:
$$\text{ENBW} = N\frac{\sum (w[n])^2}{(\sum w[n])^2} \tag{1} \label{1}$$
Where ENBW is Equivalent Noise Bandwidth (in bins), and $w[n]$ is the window samples.

The ENBW is a useful metric for window comparison and an indication of the resolution bandwidth of the window.
FURTHER DETAILS FOR THE VERY INTERESTED
The ENBW is derived from the processing gain (also called windowing loss) which is the change in signal-to-noise ratio (SNR) due to the effects of windowing (always negative compared to the rectangular window which has no loss).
The processing gain of the window is related to the ENBW as follows:
$$PG = -10\log_{10}(\text{ENBW}) \tag{2} \label{2}$$
This makes perfect sense intuitively: if the ENBW was 2 bins, then we would overestimate the total noise power by a factor of 2 (+3 dB more noise) when we sum the noise power in each bin, while the “signal” power if it only occupied one bin would not be modified relative to the noise thus resulting in a 3 dB degradation in SNR. This is detailed further below.
The processing gain is specifically due to the difference between the coherent gain that applies to the signal (when it occupies only one bin) to the non-coherent gain for the noise, from which we can derive the formula for ENBW as follows:
The coherent gain of the window, meaning the gain that would occur if all the samples of the windowed signal were in phase, would simply be the direct summation of the window, normalized by the number of samples, as follows:
$$G_c = \frac{\sum w[n]}{N} \tag{3} \label{3}$$
For example, in fred harris' classic paper http://web.mit.edu/xiphmont/Public/windows.pdf the coherent gain for the Blackman window is listed as $0.42$ which is the gain we would get as $N$ approaches $\infty$, or $20\log_{10}(0.42)= -7.54 \text{ dB}$. If we use the formula above we predict the actual coherent gain for any size $N$, such as with $N=30$ the predicted gain is 0.406 (or -7.83 dB).
>> sum(blackman(30))/30
0.406

The result of this, which provides further intuition on "coherent gain" is given in the plot below.

Likewise, non-coherent gain refers to the change in DFT output level of non coherent samples (such as white noise) due to the window function. Non coherent samples sum in power, resulting in an rms result given as:
$$G_{nc} = \sqrt{\frac{\sum w[n]^2}{N}} \tag{4} \label{4}$$
So we see that noise and signal components will change differently due to windowing and the ratio of this difference is the Windowing Loss, also called the "Processing Gain" as defined earlier and here given as:
$$PG = 20\log_{10}\bigg(\frac{G_c}{G_{nc}}\bigg) \tag{5} \label{5}$$
By equating $\ref{2}$ with $\ref{5}$ we get:
$$PG = -10\log_{10}(\text{ENBW}) = 20\log_{10}\bigg(\frac{G_c}{G_{nc}}\bigg)$$
$$= -10\log_{10}(\text{ENBW}) = 10\log_{10}\bigg(\frac{G_c}{G_{nc}}\bigg)^2$$
$$= -10\log_{10}(\text{ENBW}) = -10\log_{10}\bigg(\frac{G_{nc}}{G_c}\bigg)^2 \tag{6}$$
So
$$\text{ENBW} = \bigg(\frac{G_{nc}}{G_c}\bigg)^2 \tag{7} \label{7}$$
which by substituting $\ref{3}$ and $\ref{4}$ into $\ref{7}$ results in $\ref{1}$.
What may be confusing at first is how bins related to correlated signals can have a different power result than bins related to noise components after windowing, resulting in the change in SNR. This is explained intuitively by the ENBW: Each bin is reporting the power in its own bin plus some or all of the adjacent bins due to the spectral widening from the window. Thus in the case of white noise where all the bins are at or close to the same power level, when you sum all the bins you will over-report that actual power since power in adjacent bins is getting double-counted. In the case of a single tone (with no other tones present), it's power occupies one bin so can't get double counted (but would of course effect multiple tones due to the spectral leakage).  Without windowing (meaning using a rectangular window), the ENBW is 1 bin, so for noise the power sum of all the bins would be equal to the total power which is consistent with Parseval's theorem. This is NOT the case after windowing as explained above. We also see from this that the SNR for a waveform that itself occupies several DFT bins would also be affected differently as the waveform would thereby approach the result we get with noise, meaning reduced windowing loss since the signal and the noise would approach being affected equally.
This is also derived from $$\frac{\sigma_W+\mu_W^2}{\mu_W^2} \tag{8}$$
Where $\sigma_W$ is the variance of the window and $\mu_W$ is the mean of the window.
ENBW and PG are useful metrics when comparing window functions.
Update: I just saw this related article posted on Linked-In, I only read it over quickly but appears to be much more detailed and relevant to this post so will link here:
https://www.gaussianwaves.com/2020/09/equivalent-noise-bandwidth-enbw-of-window-functions/


Proof that spectrum analyzer measures mathematical Power Spectral Density It seems to me the power spectral density (PSD) is defined as follows. Consider a signal $X(t)$ The windowed version of this signal is

Consider a signal $X(t)$. We defined a time-windowed version of $X(t)$ as

\begin{align}
X_{\Delta t}(t) = X(t)W_{\Delta t}(t)
\end{align}

Where $W_{\Delta t}(t)$ is a time windowing function. We will choose

\begin{align}
W_{\Delta t}(t) = \frac{1}{\sqrt{\Delta t}}\theta\left(t-\frac{\Delta t}{2}\right)\theta\left(\frac{\Delta t}{2} - t\right)
\end{align}

Where $\theta(t)$ is the Heaviside function. $W_{\Delta t}(t)$ thus windows $X(t)$ to the range $-\frac{\Delta t}{2} < t < \frac{\Delta t}{2}$ and scales by $\frac{1}{\sqrt{\Delta t}}$.

The Fourier transform of a signal is defined as

\begin{align}
\tilde{X}(f) = \int e^{-i2\pi f t}X(t) dt
\end{align}

I understand the PSD to be defined as

\begin{align} \tag{1}
S_{XX}(f) = \lim_{\Delta t \rightarrow \infty} \langle |\tilde{X}_{\Delta t}(f)|^2\rangle
\end{align}

Where $\tilde{X}_{\Delta t}(f)$ is the Fourier transform of $X_{\Delta t}(t)$.

My question is how to prove than the output of a heterodyne style spectrum analyzer approximates the power spectral density.

My understanding of such a spectrum analyzer is that it takes in a signal $X(t)$ which is mixed with a local oscillator signal $f_{LO}$. This mixing brings the component of the signal at frequency $f_{LO}+f_{IF}$ down to $f_{IF}$ where the mixed signal is passed through a narrow bandpass IF filter $\tilde{F}_{IF}(f-f_{IF})$. Here $\tilde{F}_{IF}(f)$ is a filter centered about zero frequency so that $\tilde{F}_{IF}(f-f_{IF})$ is centered about $f_{IF}$. The two steps of mixing + filtering can be combined into one step of passing the signal through an effective filter $\tilde{F}_{IF}(f - f_{IF} - f_{LO}) = \tilde{F}_{IF}(f - f_0)$ where $f_0$ can be tuned by tuning the $f_{LO}$

Finally after passing through the IF filter the total power in the signal is detected.

it then mixes it to an intermediate frequency $f_{IF}$ and then passes it through a narrow bandpass filter, $\tilde{F}_{IF}(f)$ and then detects the resultant power.

Let us define

\begin{align}
\tilde{F}_{IF,f_0}(f) = \tilde{F}_{IF}(f-f_0)
\end{align}

The output of the spectrum analyzer is then

\begin{align} \tag{2}
X(t) \rightarrow (X \ast F_{IF,f_0})(t) \rightarrow |(X \ast F_{IF,f_0})(t)|^2
\end{align}

Here, because the spectrum analyzer technically measures in the time-domain, I have shown the action of the IF filter as convolution with the filter impulse response $F_{IF,f_0}(t)$.

My question is how to show $(1)$ and $(2)$ end up being the same for measured signals. It is clear that they are intuitively very similar but I'm having a hard time showing their equality mathematically.

Here are my half-baked ideas towards a solution
I know that

\begin{align}
\tilde{X}_{\Delta t}(f) = (\tilde{W}_{\Delta t} \ast \tilde{X})(f)
\end{align}

This convolution in Fourier space due to the windowing in time is very suggestively similar to the "swept filtering" which comes from demodulating the signal with the tuned LO signal. It is clear that $W_{\Delta t}(t)$ is related to $F_{IF}(f)$. However, it seems tricky because in $(1)$ the convolution happens in frequency space while in $(2)$ the convolution happens in time space.

In a spectrum analyzer the averaging happens by taking the time average of the measured power, that is

\begin{align}
X_{sig}(t) \rightarrow \frac{1}{t_{aq}}\int_{t=0}^{t_{aq}} |X_{sig}(t)|^2 dt \rightarrow \langle |X_{sig}(t)|^2\rangle
\end{align}

where in our case $X_{sig}(t)$ is the output of the IF filter. Perhaps at this point I could assume $t_{aq} \rightarrow \infty$ to apply Parsevals theorem to get something like

\begin{align}
\int |\tilde{X}_{sig}(f)|^2 df = \int |\tilde{F}_{IF,f_0}(f)\tilde{X}(f)|^2 df
\end{align}

This is starting to seem close but it's still not quite the convolution that appears in $(1)$

Can anyone give me some direction on these mathematical manipulations?

I am assuming $X(t)$ is a passband signal having finite support (non zero value for only a finite range of frequencies within spectrum analyzers range).

In your equation (1) what you have missed is $f_{LO}$ to $f_{IF}$ conversion. So $W_{\Delta t}$ should be $W_{IF}(t) = e^{j2\pi (f_{LO}-f_{IF})t}W_{\Delta t}(t)$. When you convolve $X(t)$ and $W_{IF}(t)$, in frequency domain you will get the magnitude as $|(X(f - f_{\Delta}) * W_{IF}(f - f_{\Delta})|^2$. Here $f_{\Delta}$ is the frequency shift needed to bring $X$ to IF and $W$ also includes the windowing parameter necessary to remove discontinuities around $\pm\Delta t/2$.

In (2), you are first applying a mixer-filter but then comparing (1) - a frequency domain quantity with (2) a time domain quantity. Let us say, the quantity in (2), you have applied the windowing after mixing to IF. ie, $X_{IF}(t) = (X(t)*F_{IF}(t)) \times W_{\Delta t}(t)$.

If you take fourier transform of above, you will get

$|(X(f - f_{\Delta }) F_{IF}(f - f_{\Delta})) * W_{\Delta t}(f - f_{\Delta })|^2$. But for the frequencies of interest (remember the spectrum analyzer has only finite bandwidth, so it will filter out frequencies outside it) the value of $F_{IF}(f) = 1$. So finally it will be $|X(f - f_{\Delta}) * W_{\Delta t}(f - f_{\Delta})|^2$ which is the same as we got above. Your intuition is right about (2). It is indeed convolution in time domain. But if transform it back to frequency domain, you can see it is same as (1). There are lot of things happening in the background. We have sampling and before sampling we have the Low Pass filter. All this will eliminate the unwanted frequencies outside the band of interest. Hence we can 'forget' the effect of filter $F_{IF}$


What is the theoretic difference in the ways to compute the spectrogram? I have noticed that the spectrogram for a audio signal is often calculated in different ways when the result from the DFT is transformed into real numbers; some tend to compute $|k|$ or $|k|^2$ while others use the formula $Re(k)^2 + Im(k)^2\ \ $ for each DFT bin $k$. What is the theoretic difference between these and what should be considered as the standard (if any)?

Firstly $$abs(z)^{2}=Re \left \{ z \right \}^{2}+Im \left \{z  \right \}^{2}=z*z^{^{'}}$$

These are all mathematically identical, so either one will work.

That leaves the question of whether the square should be applied or not. In nearly all cases the original signal will represent a "linear" quantity such as voltage, current, force, pressure, particle velocity, volume velocity, etc. The physical power is always related to the square (or more precise the product) of these quantities. So if you are interested in actual physical power (or energy or intensity), then you to need to apply the square. That is indeed to convention for audio signals.

[Geek Mode on]: One caveat is that simply squaring a linear quantity makes an implicit assumption about the impedance which can be misleading. For example, if you measure with a microphone reasonably far away from a loudspeaker the acoustic intensity (at this spot) is indeed proportional to the square of the measured pressure. Close up to the speaker that's not the case any more since particle velocity and pressure are not in phase. [Geek Mode off]


Basic transfer function question Why is the transfer function from node 2 to 3 $H_{23}=\frac{z^{-1}}{1+az^{-1}}$?
I don't understand the $z^{-1}$ in the numerator.


Line based Extended Kalman Filter for localization : Measurement Jacobian I have to estimate location w.r.t a given map of the environment. The state is expressed as

$ X = \begin{bmatrix} x \\ y \\ \theta \end{bmatrix}$

Given ,

$ X_{t-1} = \begin{bmatrix} x_{t-1} \\ y_{t-1} \\ \theta_{t-1} \end{bmatrix}$ and it's covariance $P_{t-1}, $ a 3X3 matrix .

I have a measurement model which detects lines representing obstacles , so , in one measurement there can be any number of lines , these lines are represented as $$ Z = \begin{bmatrix} \alpha \\ r \end{bmatrix}$$ normal form of line.

I have successfully completed all the steps until the calculation of innovation covariance i.e.

$ S = H_t\hat{P_t}H_t^T + R_t$ ,

where $H_t$ is the measurement Jacobian , $\hat{P_t}$ is the estimated motion covariance and $R_t$ is the measurement noise matrix.

When I have single line measurement , the measurement Jacobian, $H_t$ is a 2-Dimensional $2\times 3$ matrix and the measurement noise matrix , $R_t$ is $2\times 2$. But , when there are multiple line measurements , for example 8 lines , $H_t$ becomes $2\times 3 \times 8$ three-dimensional matrix and $R_t$ a $2\times 2 \times 8$ three-dimensional matrix.

Because of this change in dimension I cannot carry out the calculation of innovation covariance $S$ as $\hat{P_t}$ is a $3 \times 3$ matrix.

How can I incorporate multiple line readings and suppress the three-dimensional matrix to match innovation covariance equation ?

How are you modeling the state transitions? How are you modeling the output equation i.e. How does the state get to the measurements? Without that information, we can't tell you much about how to get your Jacobian right.

Usually, with multiple measurements, you just stack the state / observations so things go from $2 \times 3$ to $2\times 24$ in the case of 8 measurements, instead of augmenting in another dimension.


PIsarenko Harmonic Decomposition method for PSD estimation i would like to present  steps for estimation Power Spectrum density using Pisarneko Harmonic  Decomposition

i would like to present code for  step 1 part

function [vmin,sigma]=phd(x,p)

% the Pisarenko Harmonic decomposition
x=x(:);
R=Covar(x,p+1);
[v,d]=eig(R);
sigma=min(diag(d));
index=find(diag(d)==sigma);
vmin=v(:,index);
end

now let us  consider following time series

fs=100;
ts=1/fs;
t=0:ts:2.93;
x=10*sin(2*pi*10*t)+15*cos(2*pi*10.5*t)+14*sin(2*pi*10.19*t)+3*randn(size(t));

PSD using music method can be estimated  as follow

pmusic(x,6);

how can i  got the same result using pisrenko harmonic decomposition? thanks in advance

Since you have $v_{min}$, the eigenvector corresponding to the minimum eigenvalue, compute the function $\hat{P}_{PHD}(e^{j\omega})$.

The vector $e$ is given by $e = [1 \quad e^{j\omega} \quad e^{2j\omega} \ldots e^{pj\omega}]$

Make a plot of $\hat{P}_{PHD}(e^{j\omega})$ and measure the $p=3$ peaks to find the frequencies.


Behaviour of integrator at steady state I wanted to calculate response of integrator  of sinusoidal input at steady state via these two methods as mention in image  but these two methods give two different answers at steady state, so where I'm going wrong? What concept I'm missing?

Let's take a slightly more general input signal
$$x(t)=\sin(\omega_0t)u\left(t-\frac{\theta}{\omega_0}\right),\qquad\theta\in [-\pi,\pi)\tag{1}$$
It's straightforward to show that the response of an ideal integrator to the input $(1)$ is given by
$$y(t,\theta)=\frac{\cos(\theta)-\cos(\omega_0t)}{\omega_0}\tag{2}$$
The response $y(t,\theta)$ always consists of a time-dependent component $-\cos(\omega_0t)/\omega_0$ and a DC-term $\cos(\theta)/\omega_0$, which varies between $-1/\omega_0$ and $1/\omega_0$, depending on the time when the sinusoid is switched on. Consequently, no matter how large $t$ becomes, the component that all responses to switched sinusoids have in common is always the time dependent part, but not the DC-term.
The ideal response to a sinusoid extending from $-\infty$ to $\infty$ only consists of that time-dependent component that is common to all responses to switched sinusoids. You can imagine this response as an average of $(2)$ obtained by integrating over $\theta$:
$$\tilde{y}(t)=\frac{1}{2\pi}\int_{-\pi}^{\pi}y(t,\theta)d\theta\tag{3}$$
For the some discussion on the computation of the response of an ideal integrator to a sinusoidal input, take a look at this question and its answers.


Wavelet packet analysis, Sampling rate, Decomposition level I am running Wavelet packet analysis for a 1-D signal. I want to estimate dominate frequencies in the original signal.
Imagine I have following signal:

fn=[0.4 0.6 1.0];E=[0.03 0.05 0.1];A=[1.0 2.0 3.0];
wn=2*pi.*fn;a=E.*wn;wd=wn.*sqrt(1-E.^2);
y=A(1)*exp(-a(1)*t).*sin(wd(1)*t)+A(2)*exp(-a(2)*t).*...
 sin(wd(2)*t)+A(3)*exp(-a(3)*t).*sin(wd(3)*t);

so it contains these frequencies  $\big\{0.4Hz,0.6Hz, 1Hz\big\}$
My question is what is the best level of decomposition to see these frequencies?
I select $10Hz$ for sampling rate.
So my question, what sampling rate should I select for this case to identify these 3 freq?
(I've select 10Hz for sampling rate)
and what level of decomposition should I select to see all of these 3 frequencies?
Here is my little scrip for this signal with wavelet packet

df=30;N=300;                            //df:final time,N:segments
t=linspace(0,df,N); dt=t(3)-t(2);      // dt: time intervals
fs=ceil(inv(t(2)-t(1)));              // fs: sampling freq
fn=[0.4 0.6 1.0];E=[0.03 0.05 0.1];   //fn:freq, E:damping
 A=[1.0 2.0 3.0];wn=2*pi.*fn;a=E.*wn;wd=wn.*sqrt(1-E.^2);  // A:Amplitude.
y=A(1)*exp(-a(1)*t).*sin(wd(1)*t)+A(2)*exp(-a(2)*t).*...
      sin(wd(2)*t)+A(3)*exp(-a(3)*t).*sin(wd(3)*t);
level = 3;
wpt = wpdec(y,level,'coif5');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

here are result for Level 3,4,5,6 decompositions.

 as you see even in level 6 it does not show accurate frequencies. ans If I increase the level to 10 I will observe higher freq that are not exist in the actual signal.
for level 10:

"My question is what is the best level of decomposition to see these frequencies?"

You can achieve this through separation of the three components you have into different bands. Since the inner most level is a low passed version of the signal, you should try using at least 3+1 decomposition levels.

With dyadic scaling (the decomposition method typically used with wavelets) outer decomposition level covers the range of frequencies from F/2 to F .. The next level covers F/4 to F/2.. and the next level covers F/8 to F/4....etc.

Assume a signal of two sinusoids (0.4Hz and 1 Hz) with 10Hz sampling rate you get the following plot for its frequency content decomposition into 5 scales dyadic wavelet scales. by changing the sampling ratio, the location of the peaks in wavelets subbands change


Can Binary Content alone affect Noise Performance? Assume we got a transmitter and a receiver. We send a binary stream of let's say 8 bits (the stream length cannot change), using AM-modulation. We have no info on the channel/noise characteristics. Our end-system is triggered whenever a certain binary pattern arrives at the receiver (assume the triggering pattern is 01010101).

Our target is to ensure that our end-system will only be triggered whenever we send the triggering stream of 01010101 and minimize the possibility that this occurs due to noise (let's say when sending 01010100 --> receiver translates it as 01010101 and triggers the system).

Is it possible ONLY by choosing the triggering sequence appropriately to increase noise-performance? For example, to say that I choose pattern 11111111 as my triggering sequence, because it is more resistant than 01010101?

Is it possible ONLY by choosing the triggering sequence appropriately to increase noise-performance?

Yes.

The trick is to ensure that your transmitted sequences are as distant apart as possible from each other. Distant, here, is in terms of Hamming distance. This is the best you can do without more complex forward error correction codes and without any further specification regarding the nature of the noise.

What you do not mention is the nature of your message. What is it that you want to send over this link?

For example, if your application is remote control, then it is likely that you would want to transmit a small set of commands. For example, "on" and "off". This is just 2 source symbols to be encoded in 8-bit words.

A $N = 8$ bit word can describe $2^N = 2^8 = 256$ different words. The maximum (Hamming) distance between them is $8$. This is between the words $00000000$ and $11111111$. Therefore, if you map your 2 source symbols to these two words, you would have created the largest separation possible between them and it would now be difficult to send one and mistake it for the other...within reason.

This would be equivalent to a 1/8 repetition code. In fact, repetition codes that are decoded with majority voting need to be of odd length.

The more source symbols you are trying to cram into this space, the less the distance between them will have to be made and the less error-prone the whole process will be.

Alternatively, you can keep this scheme as suggested and increase your word length to 2 or more bytes in an attempt to increase this inter-codeword distance as much as possible.

Hope this helps.

EDIT:

There is nothing more that you can do to guard against noise at the level of the source message (short of using Forward Error Correction of course).

The only thing that you can do is to select words with large distances.

Now, having said this, your message might have its own statistics. This goes all the way back to the beginning.

In other words, each one of these 256 8-bit words does not have equal chances of appearing. If you were sending some text in english, for example, your most frequent letter would be "e", followed by "t", "a" and so on.

In fact, according to the second table (in the web page linked above), 17 letters of the english alphabet account for 90% of the overall frequency distribution. That is 17 out of 256 possible words have to be transmitted 90% of the times. The majority of words are made up of these letters.

What you have to do now is to select a code sequence for your trigger that is as far as possible from the codes that are used 90% of the times.

In other words, given a distribution of source symbols, the problem becomes one of "code word" mapping.

But, there is another point here that is counter-intuitive. Longer sequences have higher chances of being contaminated by noise than short sequences.

Therefore, you can take this into account during your "mapping process" by using prefix codes. Prefix codes are self-delimited. In other words, they are variable length (e.g. one sequence might be 2 bits and another 8) and no prefix sequence (that is the first N symbols of it) is part of another code. Therefore, the codes are unique and "self-delimited". If you were to see the sequence "1110100110", you know that it contains 3 symbols "111", "010", "0110". None of those is part of the beginning of another one.

Your trigger symbol could be the shortest and your message symbols could be assigned to the rest. In this way, your effective link speed will also increase because you don't have to use all 8 bits for all symbols. In some cases, you might even be able to pack 2 symbols in a byte and within reason this doubles the link speed.

Hope this helps.


IIR Filter structure implementation dynamic range I want to implement an IIR filter.  The filter is going to implement by cascading second order sections. Below there are 3 pictures that implement this filter using:

Denominator - Numerator using Tapped Delay Line Structure 2:
Numerator -Denominator using Time Delay and Accumulate Structure 3:
Denominator - Nominator using Time Delay and Accumulate Structure

(1)

(2)

(3)

a,b: coefficients

I want to examine the efficiency of the structure due to dynamic range and ease of implementation. I also want to use only one multiplier in the architecture. I think it is easier to implement the structure with less delayers but what is the dynamic range and how i can find it??

Peter's correct.  at least if fixed-point arithmetic is used.  the DF2 has poles before zeros.  that means the signal is getting boosted by the poles before the zeros (which are often sitting very close to the poles) beat the gain back down to reasonable levels.  so if the signal overflows, saturates, and distorts due to this gain, the attenuation offered by the zeros does not reverse that distortion error.  it's too late.

now the DF1 implements the zeros before the poles, and if you did not have an accumulator that has a double word width, then the DF1 suffers a sorta opposite flaw.  instead of potentially overflowing (and saturating) due to gain from the poles which is what the DF2 does, the DF1 will have attenuation first (due to the zeros) followed by the compensating gain from the poles.  if the accumulator is no wider than the signal word going in, there is roundoff error due to quantizing the attenuated signal.  then when the poles act on the attenuated and quantized signal, that quantization error is boosted by the poles.

BUT if your DSP (or whatever processor you use) has an accumulator that is twice the width of the two multiplicands, there is no roundoff error for the poles to boost.  so the DF1 avoids the overflow problems of the DF2, but if and only if there is a double-wide accumulator, does not suffer from boosting the roundoff error.


How do we use signal::filter and signal::butter in R for EEG data? I am trying to design an experiment to determine the peak amplitude of an EEG signal in response to a stimulus. Till now, our team has been using MATLAB and since we wish to go open source, we are trying to replicate the same in R.

High pass cutoff: 0.3 Hz
Low pass cutoff: 30 Hz
Number of observations: 59000
Sampling rate: 500 Sps

The code I have written is as follows:

#Trying to pass parameters as pi-radians
bf <- signal::butter(2,W=c(0.3,30)/250,type="pass",plane="s")

#test is my EEG dataframe in R
for (i in 2:9) {
  test[,i] <- signal::filter(bf,test[,i])
}

I tried following the example on stackoverflow:

How do I run a high pass or low pass filter on data points in R?

However, the amplitudes which I am getting via this is greatly different for the same set of parameters when using MATLAB's EEGLAB which uses pop_basicfilter.m.

Please let me know if you need any additional information.

I think your problem is that you want to use signal::filtfilt rather than signal::filter for the cut-offs you set to map on to the half-amplitude cut-off (otherwise, with the signal::filter function, these are the half-power cut-offs).

See my answer here for example code showing how to use the W parameter in signal::butter and double-check that your filter is working appropriately:

https://stackoverflow.com/a/49905340/9665680


How to decide what platform is best to implement real-time audio processing on? When developing products I do all of my algorithm design in Matlab. They are usually pretty basic, an IIR or FIR filter or two, a couple FFTs, etc. When it comes time to move it to an embedded environment I always have a hard time deciding what platform to run it on. I usually think of these (broad) platforms:

DSP core
FPGA
Microcontroller
ARM

What factors should I consider when trying to make this decision?

If you are keeping to "standard algorithms" like IIR, FIR, radix-2 or 4 FFT (ie stuff that fits DSP architectures well without much control flow), you can try this:

Count up how many "multiply accumulates" you need per second in all your algorithms.

< 10 million you can probably get a fast microcontroller to do the job (or even a slow one if you are <1M)
< 100M is easy DSP territory
< 1G is likely to be doable fast DSP territory
1G-10G is where the cross-over between DSP and FPGA comes
> 10G is multiple DSPs or FPGA


OpenCV Motion Analysis: Optical Flow vs Motion History OpenCV Motion Analysis Documentation lists multiple techniques for motion tracking:

Sparse Optical Flow | calcOpticalFlowPyrLK,
Dense Optical Flow | calcOpticalFlowFarneback, and
Motion History Image | calcMotionGradient, segmentMotion, calcGlobalOrientation.

I need help in understanding their pros and cons. Which ones should be preferred under what circumstances?

Context:
I am trying to track lane departure/lane changes. Any abstract/generic response would be good enough. Thanks!

Farneback is a dense optical flow algorithm. This means, it outputs the flow vectors per each pixel (which can be tracked). LK is a sparse variant, tracking only certain feature points.

The motion segmentation algorithms generally target estimation of foreground or background and doesn't allow establishing correspondences. Please refer to this post. You will find more details in my response.


Is MATLAB documentation wrong on when to use discrete versus continuous wavelet transform? I was reading some of the MATLAB documentation and I get confused when reading the part about When is Continuous Analysis More Appropriate than Discrete Analysis?

My knowledge before reading the documentation:

I always thought that a signal that is non-differentiable everywhere requires an infinite number of terms to represent it.  This is true for the Fourier basis and a wavelet basis.

The documentation:

To answer this, consider the related questions: Do you need to know all values of a continuous decomposition to reconstruct the signal exactly? Can you perform nonredundant analysis? When the energy of the signal is finite, not all values of a decomposition are needed to exactly reconstruct the original signal, provided that you are using a wavelet that satisfies some admissibility condition.

What I find confusing:

The documentation seems to imply that if the “energy of the signal is finite”, then you don’t need all the terms in the basis expansion.  I don’t think this is true, because I can think of signals with finite energy, that are not differentiable everywhere.  And if they aren’t differentiable everywhere, then they will require an infinite number of terms to represent using the Fourier basis or the wavelet basis.     A perfect example of this is the boxcar function.

Is the documentation wrong?

[EDIT] Let us start from a verbose and simple example. In $\mathbb{R}^2$, the Mercedes-Benz pattern constitute a tight frame: redundant, energy-preserving:

Each vector $x=[x_1,x_2]$ yields three three projections $[p1,p2,p3]$, from which you can recover $x$ in a number of robust ways: if the projection value are corrupted, rounded, even if one is cancelled, you can still recover $x$, at least approximately. Because of the redundancy $2\to 3$ in the representation. This is not the case with the orthogonal basis the $\{[1,0],[0,1]\}$: if you lose one projection, $x$ is almost lost (you only keep one direction).

In other words, you don't need all coefficients with a frame of vectors (or wavelets). This basic example extends to more generic frames or bases. You can find more on that within the topic of oversampled filter banks, and the paper Filter bank frame expansions with erasures. Those are used for redundant channel coding, for instance.

In the continuous setting (both the signal and the wavelet transform), due to the huge redundancy, many possible inverses exist. Notably, given an analysis wavelet, the requirements on the synthesis wavelets to verify the admissibility conditions (on a pair of wavelets) are very mild. On these topics, I would suggest two references:

Matthias Holschneider, Wavelets, an analysis tool, 1995
Bruno Torrésani, Analyse continue par ondelettes, 1995

The topic is quite involved, may involve a lot of complicated functions spaces, different for the wavelets and the signals. A wavelet analysis is given by:

\begin{equation}
W_f(b,a)\ =\ \int_{\mathbb{R}}\ f(t)\ \tfrac{1}{\sqrt{a}}\,\psi^{*}(\tfrac{t-b}{a})\ d
t\ .
\end{equation}
Provided that $\psi$ is admissible, i.e. when the
two constants $c^{\pm}_{\psi} =
2\pi \int_{0}^{+\infty}\tfrac{|\hat\psi(\pm\omega)|^{2}}{\omega} d\omega
< \infty$ are finite and equal\footnote{When $\psi$ is sufficiently
  regular, this condition reduces to a zero-average requirement, that is,
  $\int_{\mathbb{R}} \psi(t)\ d t = 0$}, that is,
$c^+_\psi=c^-_\psi=c_\psi < \infty$, the signal $f$ may be recovered
from the coefficients $W_f(b,a)$:
\begin{equation}
f(t)\ =\
\tfrac{1}{c_{\psi}}\int_{0}^{+\infty}\!\!\!\int_{\mathbb{R}}\ W_f(b,a)\
\psi_{(b,a)}(t)\ d b\,\tfrac{d
a}{a^{2}}.
\end{equation}

You can find about the same idea in Torrésani's book, page 20:

Still considering a continuous signal with finite energy, but trying to build a discretized  wavelet frame, you can sample in the wavelet space with perfect reconstruction provided the sampling is dense enough with respect to some lower bound computed from the wavelet. So effectively, you don't need all (continuous) values. But the number of samples needed could be denombrable.

Basically, you have a frame when for every unit energy signal $f$ projected onto the set of frame elements $\psi_n$:

$$A\le \sum |<f,\psi_n>| \le B\,.$$

If $A=B=1$, you recover the case of orthogonal expansions. Ingrid Daubechies (for one) showed, from the continuous formula above, that one can find some discretizations $a_n$ and $b_n$ to build frame elements

$$\psi_n = \psi(\tfrac{t-b_n}{a_n})$$

and recover $f$ from the $|<f,\psi_n>| $, sometimes with iterative formulae that converge faster if $A$ is close to $B$.

Moving on to discrete signals (where differentiability is not well-defined), you can build unions of wavelet basis, or wavelet frames like the stationary wavelet transform. In this redundant (but with a finite number of samples)  expansion, you can pick a subset of samples, that corresponds to, at minimum, a specific choice of basis vectors. Once again, you do not need all samples.

THe rationale is that if an $N$ sample signal get represented by $M> N$ coefficients, when they are well formed (frame coefficients for instance), different inverses exist, some of then discarding parts of coefficients.

For a basis (non redundant) expansion however, all samples are required in general.


Minimum number of bits required for each character Assuming that a computer is interfaced with a device that has 210 characters and
that each character is transmitted using binary words what is the minimum
number of bits needed to represent each character? To find the minimum is there any specific formula ?

the formula is

$\lceil \log_2 210 \rceil = 8$


Digital filter as physical device I'd like to find some book or books to get information about how digital filter is built depending on specifications. Like depending on price, or speed and etc. As I know it could be built by some D flip-flops and summation blocks, but there are more kind of them depending on specs, so if it's possible to get book or website about this it would be great.
Thanks in advance.

Digital Filters as dedicated pieces of hardware are quite rate these days. Most filtering just happens on general purpose CPUs such as ARM core and may be partially accelerated by specific instruction sets or co-processor (e.g. NEON). Your smart phone or your laptop do a lot of digital filtering every day.
Dedicated filter hardware is typically part of a special purpose SOC (system on a chip) that's design for a specific  application like a cable modem or ANC processor). The filter  implementations are optimized for the specific application and "configurable".


Extract a square-shaped signal which lies in a specific range I need to extract (estimate) a square-shaped signal from the main signal. This square-shaped signal should be located within a specific range, in fact all we have to do is to estimate its height. A square-shaped signal is not periodical, it's length can vary from 1/12 to 2/3 of overall length.

I've read some literature about BSS and ICA but I feel like this problem should have a very simple solution.

If you just want to estimate the height of the square shaped component you can pass the composite signal thorough a low pass filter with a cut-off frequency close to DC or less than the minimum frequency you can expect to see in the composite signal apart from the square component. The resulting signal would have most of the non-DC components filtered out and hence the energy contained in the filtered signal must be close to the that of the square pulse. After this you can do a Short time Fourier transform of composite signal with a time window less than half the minimum length of the square signal (1/24) and plot the spectrogram. Note the magnitude and variation of DC in the spectrogram. This information when corroborated with low pass filtered result should give you a fair estimate of the height and duration of the square signal.


Remove AC component from the mixed signal using an elliptic filter I have a signal like this:

$s(t) = C + a \cdot sin(2 \cdot \pi \cdot f \cdot t)$

where $C$ - arbitrary constant, $f$ - frequency of AC component (Hz) and $a$ - amplitude of AC component.

Parameters: $C = 1$, $f = 1$, and $a = 0.1$.

I want to synthesize a filter that will remove AC component with an angular frequency $\omega = 2 \cdot \pi \cdot f$ from this signal, i.e.:

$s(t) = C + a \cdot sin(2 \cdot \pi \cdot f \cdot t) \xrightarrow{After Filtering} s(t) ≈ C$

It is difficult to make this ideal, therefore, an elliptic filter was chosen as a filter, as a good frequency splitter.

After tuning it, I get the following frequency response of the filter:

The model of a computational experiment is as follows:

And the following results are obtained.

AC component has not been deleted.
Signal after the filter repeats the input. In theory, AC components with a frequency above 0.5 Hz (and in my case the frequency is 1 Hz) should be removed, and only the constant $C$ component should remain, i.e. $1$.

Where did I make a mistake?

I agree, the 1 Hz signal should be attenuated given the filter response. I suggest sweeping the input frequency to create the frequency response manually as I assume the filter was somehow scaled, possibly by not including a factor of $T$ in the mapping from the analog filter. Try a much higher frequency to see where it is nulled and then narrow in on the actual cut-off.

Further notice that the frequency response is not completely nulled at 1 Hz, but would have about 65 dB attenuation and that it also contributes significant phase distortion to the passband.

Instead of using an elliptic lowpass filter, consider using a 2nd order notch filter for eliminating AC signals. This implementation is detailed further at this post:

Transfer function of second order notch filter


Contrast stretching Equation problem How to determine $a$ and $b$ in this equation. I know how to determine $c$ and $d$ by looking at the  intensity histogram. But what about $a$ and $b$?

$$P_{out}=(P_m -c)\left( \dfrac{b-a}{d-c} \right)+a $$

I think this is the equation you are using

http://en.wikipedia.org/wiki/Normalization_%28image_processing%29

where b=newMin and a=newMax These values can also be found from your histogram. Usually you expand the new image to take up the full intensity range. The intensity range of your image is the X-axis of your histogram. If it goes from 0 to 1, then a=1, b=0 if your hsitogram is from 0 to 255 a=255,b=0

a=newMax= most likely the last value of the x-axis on your histogram
b=newMin= most likely the very first value of the x-axis on your histogram
c=Min=last  nonzero element in histogram
d=Max=first nonzero element in histogram


What is the IMAL envelope method? "Envelope methods" can be used to describe and predict the output characteristics of a highfrequency non-linear amplifier such as a travelling wave tube. Apparently one of these methods is known as IMAL.

There's several aspects to IMAL which I'd like to gain a better understanding of:

I consider an input signal, which is a superposition of sine-waves of different amplitudes, phases and frequencies. I think I know how to compute the complex envelope of this signal. Both the amplitude and the time of the complex envelope, however, are time-varying. The nonlinear characteristics of the amplifier give the output amplitude and output phase according to the input amplitude. But using matlab, I clearly cannot assign to every of my time-varying amplitudes a new amplitude - it's just too many. I assume this is where linear interpolation is needed?
Also apparently the signal is defined in frequency domain and then transformed to time domain. Why would I do that? I already know my signal is a superposition of sine-waves (which is time-domain, right?). I don't even know how to describe this signal in frequency domain.

IMAL is described in this article [1] as

  Basically, IMAL utilizes an input signal
  defined in the frequency domain. An initial random phase value is assigned to each spectral line and the spectrum is converted into the time
  domain using the inverse fast Fourier transform (IFFT). As we are interested in the intermodulation falling in the amplifier band, the band-pass input signal is then converted to a complex lowpass equivalent
  one (complex envelope) and this is distorted using the nonlinear characteristics of the microwave amplifier. A linear interpolation routine
  is used to distort the signal and it is then transformed back to the frequency domain (FFT) to obtain the desired characteristics.

[1] "Computer Simulation of Intermodulation Distortion in Traveling Wave Tube Amplifiers" by Richard G. Carter, Wolfgang Bösch, Vishnu Srivastava, and Giuliano Gatti, in: IEEE Transactions on Electron Devices, Vol. 48, No. 1, January 2001. Available online from the University of Lancaster.

Let's start with the easy one:

  Also apparently the signal is defined in frequency domain and then transformed to time domain. Why would I do that?

Because it's easy and intuitive to do that. Say you want to test with oscillations of frequencies $\mathbf{f}=\left\{f_0,\ldots,f_n\right\}\subset \mathbb R$ with (randomly chosen) phases $\mathbf{\Phi}=\left\{\phi_0,\ldots,\phi_n\right\}\subset \mathbb R$ and amplitudes $\mathbf A \subset \mathbb R,\, A > 0\, \forall A \in \mathbf A$.

So you

choose a sampling frequency $f_s > 2\min\left\{ \max (\mathbf f), f_\text{physical system bw}\right\}$,
pick a frequency resolution of $\Delta f$ that you want to achieve, and thus find $N = \left\lceil  \frac{f_s}{\Delta f}\right\rceil$
build a vector $\mathbf S\in \mathbb C^N$ such that $$
S_i = \begin{cases}
A_n e^{j\phi_n},\, n= \text{index of $f$ in $\mathbf f$}& \forall i: f=\frac{f_s}{2i} \in \mathbf f
\\
0&\text{else.}\end{cases}
$$
do an $\mathbf s := \text{DFT} (\mathbf S)$ to get your time-domain signal.

In other words, you'd think of the spectrum of the signal that you want to produce as being "empty", and then you just "put" signals at the frequencies you want them by simply setting the corresponding bins in a frequency domain window to the phase. You just "paint" the spectrum, and then IDFT it to time domain. Feels much more intuitive to me.

  I already know my signal is a superposition of sine-waves (which is time-domain, right?)

The signal is the same in both time domain and frequency domain. It's the same signal. I cannot say this too often: Time and Frequency domain are just two different bases to which you can represent a signal. The DFT and IDFT are really nothing but basis change matrices. Really.

  I think I know how to compute the complex envelope of this signal. Both the amplitude and the time of the complex envelope, however, are time-varying.

Yep, but that's just a result of overlaying multiple constant-amplitude signals. And that's what the DFT does for you.

  it's just too many

You realize that Matlab doesn't mind a couple million mathematical operations? You don't even pay per operation!

To conclude, you just generate the signal in frequency domain because it's more intuitive and much more computational efficient to do so.

Then you simply use an unspecified time-domain algorithm (of course, on the IDFT of the frequency domain signal), and then convert back into frequency domain to see the effects on the frequencies in the signal – intermodulation, generation of harmonics and possibly subharmonics etc.

I think you should get more used to frequency domain thinking – a signal in frequency domain is really 100% equivalent to a signal in time domain, but thinking about the frequencies in a signal might be much easier than thinking of the same signal as sequence of time values, especially if there's a non-zero value for every time step, but only very few non-zero frequencies.


How to know basics about convergence I apologize if the post is incorrect.

I'm a sophomore student studying Electrical engineering. As a part of an introductory course on signal and linear systems, I'm required to learn Fourier and Laplace Transformation.

But the problem is  I don't understand the concept of 'convergence' and it's really bugging me. It pops up throught the text.

So I want the suggestion to understand convergence from the beginning.

A first short shot, because I could waste hours of your time on that.

The reality (whether it exists or not) looks, sounds, appears complicated
As our senses are limited, and our brain is too slow, we need to shrink our perceptions, enough to understand most of what we need to grasp and share
So with our limitations, we need (and are only capable of) approximations of the above reality, here some signal $s(t)$
Depending on the time we have, the purpose, the extent, different kinds of approximations are required.

Convergence is the science of 'how approximations can be expected to be close enough to the reality'.

To draw a circle by hand, you only need to seize that $\pi \approx ~ 3$ (uniform 5% error). With a compas on a gridded sheet graph paper, perhaps $\pi \approx ~ 3.15$  is sufficient. For the Jet Propulsion Lab (JPL) "for interplanetary navigation, we use 3.141592653589793", see details in  How Many Decimals of Pi Do We Really Need?

Convergence is about how more complex questions are answered (see for instance how can a function can be approximated by sines or cosines) :

can a calibration function be approximated by a parabola
if I want  raise of 20 % in two years, can I have an increase of only $9.6 \%$ each of the two years (th answer is yes).

In more academic terms: suppose that there is a function $s(t)$. One has little knowledge about it, except that it wiggles. So you can choose a set of functions $\psi_\omega (t)  $ depending on a parameter $\omega$, that wiggles like exponential sines. Can we approximate $s(t)$ with a certain combination of  a finite number of $\psi_\omega (t)$? In the standard linear case, we seek a linear combination of the $\psi_\omega (t)$:

$$\sum a_\omega \psi_\omega (t)$$

The question rephrases as, for a given interval for $t$, and an increasing set $\Omega_\lambda$, $\lambda\in \Lambda$, of $\omega \in \Omega_\lambda$. The causal ideas, is: the bigger the number of terms (so the bigger the cardinal or the number of terms in  $\Lambda$). In other terms, is

$$s(t)-\sum_{a_\omega \in \Omega_\lambda} a_\omega \psi_\omega (t)$$

very close to $0$, and how, depending on $t$ and $\Omega_\lambda$. Or: how fast does

$$\Delta\left(s(t)-\sum_{a_\omega \in \Omega_\lambda} a_\omega \psi_\omega (t)\right)$$

tends to zero, where $\Delta(\cdot)$ is a measure of closeness (a distance, a norm, a divergence).


How to take the linear convolution of these two signals? How do I perform the linear convolution of the following two signals?  I am having trouble relating $x[n]$ to a series of points, like was given by $h[n]$ below.

$$x[n] = e^{j\pi n}\left\{{u[n]}-u[n-8]\right\}\quad\text{and}\quad h[n] = (-1)^{n}\left\{{u[n]}-u[n-4]\right\}$$

$x[n]$ is a finite-length sinewave of length $L=8$, and $h[n]$ is a causal filter of length $M=4$, expressed as $h[n]=\{1,-1,1,-1\}$.

The solution is:
$y[n]=\{1,-2,3,-4,4,-4,4,-4,3,-2,1\}$


SNR Signal to noise ratio of Pseudo random bit sequence (PRBS) How should I estimate the Signal to noise ratio of a 20 Gbps NRZ PRBS pattern from power spectral density--where should I take the noise floor if I can assume well the noise is AWGN?
I know roughly what the spectra should look like: https://pdfserv.maximintegrated.com/en/an/AN3455.pdf
and how an oscilloscope would take the SNR in the time domain over a distribution of samples: https://gwdata.cdn-anritsu.com/en-us/test-measurement/reffiles/Products-Solutions/11410-00919A-Eye-Diagram-AN.pdf
but how would I correctly go about using the periodogram/ power spectrum to find the SNR? Should I take the noise floor to be where out at where the spectra is flat?
I have a time domain data set of PRBS7 with noise, the data only has one non repeated segment of the PRBS7, sampled with 1ps resolution ( why the MATLAB snr() function gives it out past 250GHz) the resolution isn't so great I think because the length of this sample is 6500ps. Matlab uses a default kaiser window for filtering, I may be able to do better with different filters, but want to see if this is on the right track.

To estimate the Signal to Noise Ratio you can use the normalized correlation coefficient to the reference time-aligned sequence. I detail this approach further at this post here including the relationship between correlation coefficient and SNR. Note that "noise" in this context is any deviation of the received signal from the reference sequence such as multipath effects and band-limiting. To negate those effects modify the reference sequence to also include that or use equalization techniques on the received sequence to match the reference. In many cases these effects are second order and may be neglected.

Noise detection

If you can assume the noise is White, and you can assume the PRBS is a rectangular waveform (with a Sinc function spectrum) then the equivalent 2-sided brick-wall noise bandwidth of the Sinc with nulls at 1/T (where T is the pulse width) is 1/T. If your signal level is sufficiently in far excess of the total noise power (such as 15 dB or more) then you can estimate the signal power from total power received and compute the noise power over that equivalent bandwidth (1/T) using the estimate of the noise floor far from the signal. Pay attention to the resolution bandwidth setting on the spectrum analyzer to then determine power values / Hz (for example if you are using a 1 KHz resolution band width and a pixel displays a power level of -70 dBm, this is -100 dBm/Hz since $10Log10(1000) = 30$ dB.) For lower signal powers you can resolve the signal to noise by knowing the total power is the additive power of signal plus noise (for example if the signal is only 10 dB higher the noise, the total power measured will be 0.4 dB higher than the signal alone given $10Log(1+10^{-10/10})= 0.4$, and for 15 dB higher will be 0.14 dB from  $10Log(1+10^{-15/10}) = 0.14$)


What's the rationale behind this pulse shaping function? Section 13.2 of the HD Radio standard for AM defines a pulse shaping function that appears to be some strange variant on a raised cosine filter:

More specifically, H(ξ) is a raised cosine filter with total width (1+α)T, and W(ξ) appears to be a "stretched" version. Why do they do this?

The pulse envelope is convolved with a Gaussian function to soften the transition between the symbols. If the standard pulse envelope were used without any softening, it would be audible to an AM detector (it runs at ~172 Hz).

Source: The IBOC Handbook: Understanding HD Radio (TM) Technology by David P. Maxson


MATLAB, How to filter a discrete signal? I have a 1-by-10000 matrix of double`s stored in a file named "fecg.mat". The matrix represents the magnitude of a recorded FECG signal.

I've plotted it against time (from 0 to 9999):

For removing the baseline I wonder if I can use a high-pass filter. How do I design a proper filter?

P.S
Signal processing isn't my field of study. I have no idea how to filter a discrete time-domain signal.

Easiest way to remove baseline is to remove the average:

filtered = original - mean(original);

Actually, the average is the first coefficient of the Fourier transform, so it is actually a very sharp filtering: you are eliminating the DC coefficient.

If you want more classic filtering, then check functions like butter and friends, which will synthesize an IIR filter, then use filter to filter out your signal.

Matlab also includes a filter design tool.


Spatial multiplexing MIMO, antenna critical spacing, antenna isolation How does antenna isolation affect single user spatial multiplexing MIMO performance in an antenna array? Text books talk about a critical antenna spacing of 1/2 of the wavelength  provided that the environment is  a rich scattering one . Would that be the sole criterion? or would the antenna-antenna isolation rather be the criterion to observe?

What you ultimately want is that the channel coefficients are uncorrelated. MIMO (and diversity techniques in general) depend on the fact that it's unlikely that many different channels will all be bad at the same time. When the channels are correlated, the performance of your receiver will go down dramatically, because correlation means that when one channel is bad, the others are likely to be bad too.
Regarding spacing: if the antennas are too close together, then they see essentially the same channel. So, some spacing is needed; the more distance between antennas, the better. The $\lambda/2$ antenna spacing is a rule of thumb that may or may not work in practice.
Regarding isolation: It turns out that when antennas are close together, they become "mutually coupled". The voltage produced by each antenna depends not only on the received signal, but on the voltages produced by the other antennas. Mutual coupling is reduced with increased spacing, or it can be attenuated using "decoupling" circuits.
Both of these effects result in an increase in channel correlation, and a decrease in performance. Most textbooks mention only spacing, and not coupling, even though this effect has been known in the RF community for a long time. See for example "Mutual coupling effects on the capacity of multielement antenna systems", by T. Svantesson et al.
To answer your question: in an actual system implementation, both spacing and coupling must be taken into account. However, keep in mind that the ultimate objective is to have decorrelated channels; antenna spacing and decoupling are just two techniques to achieve that.


Adding Noise in Frequency Domain vs Time Domain I am trying to analyse the effects of SNR level on my data.
I have a set of signals in frequency domain with f_min = 0.5GHz and f_max = 10.5GHz centered at f_c = 5.5GHz. These are simulated signals and are noisless (ideal data). Using the awgn function of MATLAB, I add noise of various SNR levels (SNR = 10, 20, 30, 50 and 60dB) in frequency domain. I now need to add the same level of SNR to their time-domain counterparts.

I came across something which I didn't quite understand. It goes as follows :
SNR_freq is computed from the data over the whole frequency,
while in practice, SNR_time matched to pulse data and is computed from the peak of the signal magnitude (A_max) as shown below :
SNR_freq = 10*log10(E/sigma^2)
SNR_time = 20*log10(A_max/sigma)
where sigma is the variance of the signal.

what exactly is the relationship between SNR_time and SNR_freq?

I usually like to add noise with my own code.
Adding Whit Noise is really easy task given the signal as you have.

All you need is to calculate your signal second moment at the frequency and add noise to the frequency bins such that the second moment of the noise creates your desired SNR.

Since the DFT is unitary transform, adding white noise at frequency domain is equivalent to adding noise at time domain.

There 2 options here, are you looking for adding White Noise (Namely spread all over your frequency) or Band Limited (Colored) noise?

For white noise just add Noise to each bin in the frequency domain such that the noise Second Moment (The Means squared multiplied by the number of bins) is as your desired SNR.

For colored noise do the same but only for the relevant frequency bins (Where the second moment is calculated only on the number of "Active" bins).


Causality, Region of Convergence and which side is the system on I want to try to summarize some relationships between z transform, ROC and causality. But I am getting confused between left/right sided sequences, the ROC and causal.

If you need stability then the ROC must contain the unit circle. I can understand that.
If you need a causal system then the ROC must contain infinity and the system function will be a right-sided sequence. I am fine with  that too.
If you need an anticausal system then the ROC must contain the origin and the system function will be a left-sided sequence. Yes that's fine.
If you need both, stability and causality, all the poles of the system function must be inside the unit circle.

But #4, doesn't it violate number 2 because we don't include infinity?

Got it from here: http://en.wikipedia.org/wiki/Z-transform#Example_2_.28causal_ROC.29

Thanks

If you combine the two requirements (stability => the unit circle must be inside the ROC), and causality (the ROC is given by $|z|>r$ with some real-valued radius $r$), then it follows that $r<1$ must be satisfied. This means that the ROC is outside a circle with radius $r<1$, and, consequently, all poles must be inside this circle with radius $r<1$. So all poles are inside the unit circle (because, by definition, there are no poles in the ROC).


Non human formant frequencies Some animals are specially trained to repeat human phrases, others just make sounds that people interpret as speech.

In this video a rooster seems to say "Cristo Viene" (Christ is coming).

Taking the spectrogram of that rooster to analyze the formant position of the vowels I can't find any matches (I was using a formant position UCLA's chart from   here ), because the pitch of the animal is obviously different to formant position of human male or females.

Does anybody knows the way to find the formants position at different pitches? (I put the spectrogram of the rooster above)

I think you have a misunderstanding about formants, pitch and partial tones. The spectrogram shows the position of the partials and you can also estimate their fundamental frequency, or pitch, trajectory to be the position of the lowest partial. The formant structure is not readily recognizable from this representation.

Formants are the resonances of the vocal tract but these resonance frequencies are not necessarily excited because they don't coincide with the partial frequencies that are at integer multiple positions of the fundamental frequency. Rather, the formants shape the amplitude of the exciting signal, in this case best described by the partials. You can think of the formants as the peaks of the (smooth) envelope that you can draw on a frequency-amplitude diagram showing the amplitudes of the partials.


How to compare 2 matrices of MFCCs, keeping the time information I have a corpus of several thousand sounds. Now I pick a sound sample S, and I want to find in my corpus the sound S' that is the most similar to S.
The application is for music, for a composition that I am working on :) I found a lot of questions which are all related to speech recognition. The problem here is different as rhythm is the most important info for me, even before timbre. If S' has a different length than S, it must be somewhat of stretched version of it.

The samples are between 0.5 seconds to 2 seconds (approximately). I have tried to compute MFCCs, apply DTW, but I don't get convincing results. I still believe MFCCs are the way to go, however I started to think DTW is not suited for my purpose, because it removes the rhythm information. For example (and if I understood well), DTW would recognize that the sounds "Puuuum pum tchak" and "Pum pum tchaaaak" are the same. But obviously, when it comes to music, they are not, as the timing is crucial.

So, how to compare 2 matrices of MFCCs, keeping the time information? Or is my approach completely wrong?

I was currently looking into procrustes analysis, as I have read that it compares 2 shapes, removing scaling, ... and this sounded good to me.

EDIT

With DTW I might have :

distance("Pum pum tchak", "Pum pum tchaaaak") < distance("Pum pum tchak", "Puuum puuum tchaaak")

in other words, the time wraping doesn't care about the fact that "Puuum puuum tchaaak" is just exactly the same as "Pum pum tchak" but with slower tempo.

You are misunderstanding how DTW works. You can actually adjust the frame deletion/insertion costs to adjust how much temporal mismatch you can tolerate.

We would have: DTW("Puuuum pum tchak", "Puuum pum tchak") < DTW("Puuuum pum tchak", "Pum pum tchk") < DTW("Puuuum pum tchak", "Purr purr purr").

That is to say, if there is a sample with a perfectly good timing, it'll always yield a higher DTW score than a sample with the wrong timing. The trade-off between similar timbre and similar timing can be adjusted by tweaking the DTW alignment costs. One useful thing with DTW is that it finds the optimal alignment between the "query" and the found sample - so you can use the path it founds to timestretch/phase-vocode the result and align it perfectly with the query.

It is not clear what kind of sensitivity on timing your application has; which features matter the most (overall timbre? onomatopoeia which is pronounced?); and what kind of processing you could do to the found clips to align them even better with the query. If rhythm/timing is absolutely essential, you might want to process each sample by an onset detector ; and restrict your search to the set of clips which have the right number of onsets, and the correct inter-onset times.

If you have absolutely no tolerance on differences in rhythm, then you could very well compute the average of pairwise distances between the two clips (padded with zeros if their length does not match). That would be a bit silly, but your application might want that.

Finally, if both your query and your samples can accept the same domain-specific representation (say a beatboxed rhythm and a drum loop can both be transcribed into "bass drum" and "snare drum" events) - it might help to transcribe all clips and to match the transcriptions instead of the audio signals.


Upsampling by zero padding or duplicating sequence I need to upsample a 2 MHz signal by a factor of 3000 to convert the OQPSK modulated baseband signal to passband. When I upsample by duplicating the elements in the sequence, I seem to get less error compared to padding the sequence with zeros. The frequency spectrum is also better when I duplicate elements whereas with zero padding, the spectrum shows nothing. I am now confused on which way is the right method to upsample a sequence and why.

Matlab’s ‘upsample()’ command does not “pad” a sequence with zero-valued samples.  The ‘upsample()’ command “stuffs” a sequence with zero-valued samples.  “Zero padding” and “zero stuffing” are two different operations.  “Zero padding” means appending a sequential string (a sequence) of zero-valued samples to the beginning or end of a sequence.

I believe what you want to do is: ‘Upsample()” your time-domain sequence to the desired final sample rate that's compatible with your 2.4 GHz carrier sequence.  Next, lowpass filter the upsampled sequence.  Then multiply the filtered sequence by your high-frequency carrier sequence.


homography matrix problem i have mini project to do about homography , but i am really a beginner in matlab . i have calculated the homography matrix easily by taking a set pixels from both images but the problem is how to apply this matrix to all pixels of the original image (i think it's difficult to treat each pixel individually ) . here is the original image :

i want to apply the homography to this image ( i have already extracted the homography matrix ) to get a frontal view of the cover.

so with my weak knowledge of matlab , i applied the folowing operations:

i get 4 points from the corners of the original image and the corresponding points in the output image

i have estimated the homography :

0.9638 -0.0960 52.5754;
0.2449 1.3808 -17.0081;
-0.0001 0.0013 1.0000

then i have applied it on the image using this poor matlab code:

function [Ired] = affect(I,P)
I1=im2double(I);
for i=1:308
i
for j=1:288
v1 = [j;i;1];
v2 = P*v1;
v3 =v2/v2(3,1);
if (v3(1,1) <0 || v3(2,1)<0 || v3(3,1)<0)
   continue;
else
Ired(round(1+v3(2,1)),round(1+v3(1,1)))=I1(i,j);
end
end

end

and the result is the following:

- but the problem with this code is it takes a lot of time (sometimes 10 minutes) to perform the operation on the image especially if the size of the image is high so is there a way to implement this code without loops to make it faster? .

- i get these black dots all over the image , i can do a dilation but the image will be distorted , do you have an idea how to eliminate this dots without distorting the image ?.

Matrix/vector operations are really strongly optimized in Matlab/Octave. Use them whenever you can.

In your case, instead of multiplying the 3x3 homography matrix with a 3x1 vector for NxM times you can easily modify your calculation to do a multiplication with a 3x3 matrix and a 3xM matrix for N times.

Check out this Octave code:

warning ("off", "Octave:broadcast");
w = 1280;
h = 720;
H = rand(3);

tic();
for i=1:h
    for j=1:w
        H*[i;j;1]/(H(3,:)*[i;j;1]);
    end
end;
toc();

tic();
for i=1:h
    v=[ones(1,w)*i;1:w;ones(1,w)];
    (H*v)./(H(3,:)*v);
end;
toc();

It gives me the following result:

Elapsed time is 25.5651 seconds.
Elapsed time is 0.273893 seconds.

Maybe you can go further by generating a 3x(N*M) matrix with all the image coordinates.

There is a simple solution to the 'black dots' problem. Right now you have the homography $H$ between the original image and the cover image and you generate the points by looping through the pixels of the original image $(x,y)$ to color the points $(x',y')$ in the destination image.

Let's have the homography $H^{-1}$, loop through all the pixels in the destination image and calculate the pixel address of the original image you want to sample the color from.


Fourier transform is an isomorphism...but we don’t get when each frequency appears? Statistician here who wants to get some DSP knowledge for time series analysis.

I’ve known for years that if we hit a function with a Fourier transform, we have an inverse Fourier transform that will recover the original function. However, doesn’t the interpretation of the Fourier transform in the frequency domain lack the time component? In other words, we can say that 100Hz appears in the signal with some intensity, but we can’t say if that appears at the beginning of the signal or the end. Those are very different signals to me, yet I’m supposed to be able to recover each by applying the inverse Fourier transform?

There seems to be an inconsistency: we either lose the time information and can’t invert, or we retain the time information and can invert.

What is the resolution to this apparent inconsistency?

(I have a hunch that it has to do with the imaginary part of the Fourier transform, though I’m not sure how.)

It's true that taking the Fourier transform will leave you without any (visible) information on time and vice versa, but of course you don't lose any information, you just represent it in a way such that in one domain you only see time information, and in the other you only see frequency information.

Take as an example the Fourier transform of a time-inverted (real-valued) function $x(-t)$: its Fourier transform has the same magnitude as the Fourier transform of the original function $x(t)$. The difference between these two Fourier transforms lies exclusively in the phase. So you're right to assume that timing information is encoded in the phase of the Fourier transform.

There is no time localization in the Fourier transform because its basis functions are complex exponentials extending from $-\infty$ to $\infty$. There are other transforms that will give you a certain degree of time and frequency localization, the most well-known of which is probably the Short-time Fourier transform. Also take a look at this related question and its answers.


How to analyse anti-causal discrete transfer function using matlab? Consider a discrete transfer function that represents an anti causal filter such as a derivative filter:
$$H(z) = (-z^{-2} -2z^{-1} +2z +z^2) (1/8T)$$

Where T is the sampling period.
Normally in MATLAB, I use fvtool and enter coefficients but these coefficients are relevant to negative exponentials only, i.e. for a causal filter.

I want to view the magnitude response, phase response, pole zero map of an anti-causal filter in MATLAB. How is that possible?

The transfer function
$$H(z)=-z^{-2} -2z^{-1} +2z +z^2$$
can be written as
$$\begin{align}
H(z)&=z^2z^{-2}\left(-z^{-2} -2z^{-1} +2z +z^2\right)\\
&=z^2\left(-z^{-4} -2z^{-3} +2z^{-1} +1\right)\\
&=H_1(z)H_2(z)
\end{align}$$
 It can now be seen more easily that

$H_1(z)$ has two zeros at $z=0$ (and two poles at $\pm\infty$) .
$H_2(z)$ has four poles at $z=0$ and four zeros: one at $z=1$ and three at $z=-1$. What can be said about pole-zero map of $H(z)$ overall?
Since $|H_1(e^{j\omega})|=1$, we cnclude that $|H(e^{j\omega})|=|H_2(e^{j\omega})|$
Since $\angle H_1(e^{j\omega})=2\omega$, we cnclude that $\angle H(e^{j\omega})=\angle H_2(e^{j\omega})+2\omega$


Reconstruction of a signal using 1D discrete wavelet There is a signal of $50\textrm{ Hz}$ and $120\textrm{ Hz}$ corrupted with noise. The sampling rate is $1000\textrm{ Hz}$. Here I used a 3-level DWT to extract this two components of the signal respectively. The figure is the power density spectrum of signal reconstructed from the detailed coefficient. My question is, why there is a unknown frequency component ($129.9\textrm{ Hz}$) in it?

clear all
Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1024;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')
f = Fs*(0:(L/2))/L;
h=boxcar(L);
P1=psd(X,L,Fs,h);
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
ylim([0 300])

[c,l]=wavedec(X,3,'bior4.4');
for i=1:3
d{i}=wrcoef('d',c,l,'bior4.4',i);
a{i}=wrcoef('a',c,l,'bior4.4',i);
Pa{i}=psd(a{i},L,Fs,h);
Pd{i}=psd(d{i},L,Fs,h);
end

figure
plot(f,Pa{3})
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P(f)|')

figure
plot(f,Pd{3})
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P(f)|')

This sounds like a typical effect due to downsampling after a convolution with a non perfect lowpass, when one of the frequency in the approximation becomes close to the cutoff $f_c$ at some scale, of the form $f_c = f_s/2^j$ for some $j$. Notice that the sum of $120+129.9 \approx 2*125$. Set the second frequency to $124$, you can see a second peak at $126$ Hz, similarly.

So this is aliasing, resulting from a two-fold subsampling after a FIR lowpass, whose transition band is not sharp enough. This can be avoided be sharper wavelets (like $M$-band wavelets), or resorting to shift-invariant or stationary wavelets (swt.m instead of dwt.m). The following real-life example come from seismic data processing: the data is corrupted by a $60$ Hz powerline sine (US standard), and acquired at $4$ ms. The top graph is a seismic data spectrum, after the data has been processed by the discrete wavelet, and shows a second peak at $65$ Hz ($60+65=125$), while the second one is processed by a SI wavelet.

Your question revived my memories from 2002, Efficient Coherent Noise Filtering, An application of shift-invariant wavelet denoising (source of the picture).


Testing for Linearity and Shift-Invariance I started school a week and a half ago and we have gotten already a big chunk of material and I wanted to check if what I am currently doing is correct. Here below are the problems to work out that help us for the exam in the future. (This class is signal processing)

for part (b)

Your demonstration of shift-invariance is not correct.
You time-shift the input signal by $n_0$, which is correct. But you also time-shift the sinusoidal, and this only happens if you were time-shifting the output.
The description of the system contains a signal that explicitly depends on time ($\sin(2\pi/7 n \cdots)$). In general this is a red-flag for non time-invariance.


morphological processing : Hit or Miss shape detection I am currently looking at shape detection for Hit or Miss. I am having a bit of trouble with the concept of the window that they have you to detect the shape. Up till now reading about hit or miss they just give an example of a target set of pixels and a target image 1 and the step for this are as follows:

The first step in Hit-or-miss translation is the erosion of the image by the target configuration (shape), thus producing the hits. The logical operation is given by $A\ominus B$ where A is the image and B is the target shape to be found.
The second step, we take the complement of the image $\bar{A}$ and the target shape $\bar{B}$. The we erode the compliment to produce the misses with in the image. The logical operation for this is given by $\bar{A}\ominus\bar{B}$.
The final step take the intersection of the two eroded images and any surviving pixels give the location of the target shape in the original image A. The logical operation for this final step is given by $(A \ominus B) \cap (\bar{A} \ominus \bar{B}$).
Any surviving pixels give the location of the target shape in the original image A

However when I looked further into this method for shape detection I came across the following website:
(http://oramos.epizy.com/teaching/16/dip/lectures/07_Morphological_Image_Processing_1.pdf?i=1)
and on slide 34 to 35 it has the following for shape detection

in there example for the shape detection they are using another structuring factor which they have named window, which it what I am having and issue with. with why have they done this? Is this a generally way when using hit or miss for shape detection?

They are using the same method. The second resource is just showing the windowing function. Look at slide 12, the windowing function is shown here.

Erosion implies you use a windowed kernel, and slide it around the entire image (as seen in this gif), so the first source doesn't mention this step. The 2nd source is explicitly stating you must pad the search space, and then slide the kernel around it.Hope that makes sense


What does it mean for the DFT phase to be relative to a cosine wave? The following paragraph from Understanding Digital Signal Processing got me puzzled:

  The answer is: The DFT phase at the frequency $mf_s/N$ is relative to a cosine
  wave at that same frequency of $mf_s/N$ Hz where $m = 1, 2, 3, ..., N−1$. For example, the phase of $X(1)$ is −90 degrees, so the input sinusoid whose frequency is $1\cdot f_s/N = 1000$ Hz was a cosine wave having an initial phase
  shift of −90 degrees.

I have no idea how the DFT phase is related to a cosine wave in the way just mentioned above. Any hints would be greatly appreciated.

A cosine wave is an even function around the origin (X axis = 0).  The basis vectors of the real components of a DFT are all cosine waves that start at 1 at the first sample.  Imaginary components of DFT basis vectors are all sine waves that are odd functions around the origin 0, and start at zero at the first sample in a DFT window.

However, for strictly real input, you can add a phase offset p to a cosine wave function to make cosine(x+p) look like any ratio mix of real and imaginary basis vectors (cosine waves and sine waves) in the DFT result, and that phase offset is what this book is calling the DFT phase.

For even length FFTs of non-fictional data, I prefer to do an fftshift to move the phase reference to the middle of the FFT window, instead of leaving the phase reference at the potential circular discontinuity at the first sample of the DFT window.


Help with Determinants Okay so i understand the first det, which is why i got the correct solution. However for the second det, im not sure how to solve it since because of the different coefficents (0,2,4) and i tried of different solutions. If anyone can give me a clear way to solve these it would be appreciated. Thank You.

If you open the determinants along the middle column, we get for the first determinant:

$ -(bf-ec) + (af-cd)  -(ae-bd) = -3$ $\hspace{1cm}$ $ [1]$

And for the second determinant, we get

$-(bf-ec) + 2(af-cd)  -3(ae-bd) = -2 $ $\hspace{1cm}$ $ [2]$

The required determinant is

$  2(af-cd)  -4(ae-bd) $ $\hspace{1cm}$

Which is just $2*([2] - [1])$, which boils down to $2$.


Cleaning telephone audio I have telephone audio at 16kHz and I am trying to use it with a speech transcription engine that converts spoken letters to text. However, the engine has some difficulty transcribing certain letters such as "p" and "t" due to the low quality of the audio.

This is being used in a real time application so although I would like the user to speak clearly, into the microphone in a quiet environment, this won't often be the case. The transcription engine has its own background noise removal however it isn't trained on telephone audio.

In light of my ignorance, I was hoping someone could make some suggestions on how to artificially "improve" the quality of the audio. I know this is vague but I would like to try and make the speech in the audio more discernable and experiment with the transcription engine.

Can I try and remove the static sound?
Is it possible to "amplify signals" that are unique to a certain letter?

Thanks

The transcription engine has its own background noise removal however it isn't trained on telephone audio.

There are many engines which are trained on bad quality audio and which perform just right.

It is pretty hopeless to repair the quality when it is already lost, it is better to find an accurate engine.


FFT of signal with starting value zeros I have a signal,

fs = 100;                                % Sample frequency (Hz)
t = 0:1/fs:10-1/fs;                      % 10 sec sample
x = (1.3)*sin(2*pi*15*t)                 % 15 Hz component
+ (1.7)*sin(2*pi*40*(t-2)) ;             % 40 Hz component

when I take fft of x then its amplitude and phase spectrum are:

And now if deliberately I will put the leading 10 values to zeros

x(1:10) = 0;

then as per my understanding the frequency of the complete data set is not changed so the fft would not be changed. Which is supported by the amplitude fft response.

But the phase plot is showing some strange result.

Can some help me why the response is changed. Although they both are following the same equation. Also If I have a equation cos(2*pift) and my t is from 3:Tp and before time =3sec the response is zero, then will it not have the same response as t:0:Tp.

I expand and modify your code:

N = 100;  %number of sample
fs = 100; % Sample frequency (Hz)
t = 0:1/fs:(N-1)/fs;
x = sin(2*pi*10*t);  % 10 Hz component
%x(:)=1;
f = 0:fs/N:(N-1)*fs/N;

subplot(2,2,1);
plot(t,x);
title('full signal'); grid on;
subplot(2,2,3);
plot(f(1:N/2),abs(fft(x)(1:N/2)));
title('spectrum'); grid on;

x1=x;     %x1 - zero pad signal
x1(1:ceil(N/3))=0;
x1(ceil(2*N/3):N)=0;

subplot(2,2,2);
plot(t,x1,'b',t,x1,'*');
title('zero pad signal'); grid on;
subplot(2,2,4);
plot(f(1:ceil(N/2)),abs(fft(x1)(1:ceil(N/2))));
title('spectrum'); grid on;

You can see signals before and after zero padding are quite different. And the spectrums are quite different too. There is strong part of signal x1 with frequency 10 Hz, and you can see strong peak on spectrum.

You can uncomment line with x(:)=1; First full signal will be DC (constant), second (zero padding) - rectangle. Quite different signals and quite different spectrum.

UPDATE
About phase spectrum. Update to code:

F=fft(x);
angle(F(11))  % phase = -1.5708
F1=fft(x1);
angle(F1(11)) % phase =  -1.6018

Phase for x makes sense only for bin with non zero amplitude. Non zero amplitude numbers are 11 and 91. Phase for any other index are results of calculations errors of FFT function.
As x1 has  big part of x the phase spectrum for this bins has approximately the same value (as you can see).

Hint. You can present x1 as product of $sin$ function and rectangle function. So spectrum of x1 is a convolution of spectrum $sin$ and rectangle.


Implementation of IIR Filter in dsPIC 33EP - Q15 I am implementing an IIR controller in dsPIC 33EP (16-bit, fixed-point, two's complement, wrap-around, 40-bit accumulator, 32-bit multiply, 12-bit ADC) and I am using cascaded, Direct Form I second order sections (ordering up).

I am using Q15 (fractionals), but my coefficients are within the range +/-1. The question is: If a coefficient is exactly +1, will it fit in Q15?

The confusion is due to the fact that Q15 allows values from -1 to 0.99999999906 only. So I am afraid it will be quantized to 0 or -1 in 2's complement wrap-around.

I would very much appreciate any help.

Ben

There's an old trick you can use:

$$1\cdot x=x$$

So you don't need to multiply with a factor that equals $1$.


Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers I have been working on the topic of camera pose estimation for augmented reality and visual tracking applications for a while and I think that although there is a lot of detailed information on the task, there are still a lot of confussions and missunderstandings.

I think next questions deserve a detailed step by step answer.

What are camera intrinsics?
What are camera extrinsics?
How do I compute homography from a planar marker?
If I have homography how can I get the camera pose?

It is important to understand that the only problem here is to obtain the extrinsic parameters. Camera intrinsics can be measured off-line and there are lots of applications for that purpose.
What are camera intrinsics?
Camera intrinsic parameters is usually called the camera calibration matrix, $K$. We can write
$$K = \begin{bmatrix}\alpha_u&s&u_0\\0&\alpha_v&v_0\\0&0&1\end{bmatrix}$$
where

$\alpha_u$ and $\alpha_v$  are the scale factor in the $u$ and $v$ coordinate directions, and are proportional to the focal length $f$ of the camera: $\alpha_u = k_u f$ and $\alpha_v = k_v f$. $k_u$ and $k_v$ are the number of pixels per unit distance in $u$ and $v$ directions.

$c=[u_0,v_0]^T$ is called the principal point, usually the coordinates of the image center.

$s$ is the skew, only non-zero if $u$ and $v$ are non-perpendicular.

A camera is calibrated when intrinsics are known. This can be done easily so it is not consider a goal in computer-vision, but an off-line trivial step.

Some links:
ftp://svr-ftp.eng.cam.ac.uk/pub/reports/mendonca_self-calibration.pdf

What are camera extrinsics?
Camera extrinsics or External Parameters $[R|t]$ is a $3\times4$ matrix  that corresponds to the euclidean transformation from a world coordinate system to the camera coordinate system. $R$  represents a $3\times3$ rotation matrix and $t$ a translation.
Computer-vision applications focus on estimating this matrix.
$$[R|t] = \begin{bmatrix} R_{11}&R_{12}&R_{13}&T_x\\R_{21}&R_{22}&R_{23}&T_y\\R_{31}&R_{32}&R_{33}&T_z \end{bmatrix}$$
How do I compute homography from a planar marker?
Homography is an homogeneaous $3\times3$ matrix that relates a 3D plane and its image projection. If we have a plane $Z=0$ the homography $H$ that maps a point $M=(X,Y,0)^T$ on to this plane and its corresponding 2D point $m$ under the projection $P=K[R|t]$ is
$$\tilde m = K \begin{bmatrix} R^1 & R^2 & R^3 & t \end{bmatrix} \begin{bmatrix} X \\ Y \\ 0 \\ 1 \end{bmatrix}$$
$$= K \begin{bmatrix}R^1&R^2&t\end{bmatrix} \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix}$$
$$H = K \begin{bmatrix}R^1 & R^2 & t \end{bmatrix}$$
In order to compute homography we need point pairs world-camera. If we have a planar marker, we can process an image of it to extract features and then detect those features in the scene to obtain matches.
We just need 4 pairs to compute homography using Direct Linear Transform.
If I have homography how can I get the camera pose?
The homography $H$ and the camera pose $K[R|t]$ contain the same information and it is easy to pass from one to another. The last column of both is the translation vector. Column one $H^1$ and two $H^2$ of homography are also column one $R^1$ and two $R^2$ of camera pose matrix. It is only left column three $R^3$ of $[R|t]$, and as it has to be orthogonal it can be computed as the crossproduct of columns one and two:
$$R^3 = R^1 \otimes R^2$$
Due to redundancy it is necessary to normalize $[R|t]$ dividing by, for example, element [3,4] of the matrix.


If $X$ and $Y$ are zero mean independent Gaussian random variables with different variances, what is the density of $\sqrt{X^2+Y^2}$ Let $X\sim\mathcal{N}(0,\sigma_X^2)$ and $Y\sim\mathcal{N}(0,\sigma_Y^2)$ be independent Gaussian random variables. What will be PDF of $Z=\sqrt{X^2+Y^2}$ and $W=\arctan{\left(\frac{Y}{X}\right)}$. Will they still be Rayleigh and uniform distributions respectively?

According to your description, $x = Z \cdot \cos(W),\ y = Z \cdot \sin(W)$. Follow this answer for the derivation of joint PDF of $(Z, W)$ : Complex Gaussian Magnitude and Phase Joint PDF Derivation

You will reach the following expression after Methid of transformation:
$$f_{Z,W}(z,w) = |\mathbf J|.f_{X,Y}(z \cdot \cos(w), z \cdot \sin(w)), \ where \ \mathbf J \ is \ Jacobian$$
Since the computation of $\mathbf J$ will not change hence:
$$\mathbf J= z$$

Now, when you put the joint PDF of 2 independent Gaussian random varaibles $X, Y$ in the above expression you get the following :
$$f_{Z,W}(z,w) = |\mathbf J|\cdot f_{X,Y}(z\cdot \cos(w), z\cdot \sin(w)) $$ $$= z \cdot \frac{1}{\sqrt{2\pi \sigma_x^2}}\exp(-\frac{1}{2\sigma_x^2}(z\cdot \cos(w))^2)\cdot \frac{1}{\sqrt{2\pi \sigma_y^2}}\exp(-\frac{1}{2\sigma_y^2}(z\cdot \sin(w))^2)$$

Combine this and find the marginal PDFs of $Z$ and $W$ from the above expression. You will find that integrating from $-\pi$ to $\pi$ w.r.t. $w$ will not give you Rayleigh Distribution anymore. And, similarly, integrating from $0$ to $\infty$ w.r.t. $z$ will not give you $\frac{1}{2\pi}$, and hence not Uniform Distribution anymore.


Classifying this signal in dimension and number of channels The question is:

"Weight and height measurements of a child taken every month"

Classify this signal as:

one/multidimensional
single/multichannel

I'm not understanding how this can be described as a signal and hence how can it be classified.

I don't feel like it's a very well-posed question, as multi-channel signals are sometimes represented as vector-valued, which would make them multidimensional. However, I assume the answer that is being fished for is:

This is a one-dimensional, two-channel setup. The two channels are the child's height and weight. Each signal is one-dimensional, measured across time. The sample rate of each signal is $1\text{ month}^{-1}$. You could state them mathematically as:

$$
h[n] \equiv \text{ height of child at month } n
$$
$$
w[n] \equiv \text{ weight of child at month } n
$$