The cross pattern is typically a border effect, due to the periodicity induced by the standard implementation and hypotheses behind the Fast Fourier transform, when the image lacks periodicity from the right to the left, and the bottom to the top. In other words: if two opposite borders lacks continuity in values (when glued together), artifacts show.

The easiest way out is to multiply your image by an apodizing window, that goes smoothly to $0$ on the border (imagine for instance a 2D Hamming window).
A mirror extension is a second option.
More involved, you can try the Periodic plus smooth image decomposition, that has been used in random texture synthesis. Some Python versions have been offered:

Periodic plus smooth image decomposition, as described by Moisan
Python implementation of Moisan's periodic-plus-smooth decomposition of an image, with an improvement described in a series of blog posts, like: On the periodic-plus-smooth decomposition of an image, part 7: improved implementation of Moisan’s algorithm.


Averaging by N but keeping the same sample rate is equivalent to filtering input "y" with an FIR filter "B". Where B = 1/N * ones(1,N);

yFiltered = filter(B,1,y);

Now you simply need to downsample the data by N. I'm not sure exactly how the engineer implemented the original algorith, but I would do this.

yDownsampled = yFiltered(N:N:end);

Basically start at index "N" and keep 1 out of N samples until the end.

Your starting index can be anything between 1 and N, you can toy around with it.

You'd need to do the same with your time vector.


When can we write Heisenberg uncertainty Principle as a equality?

It is important to define the time and frequency widths $\Delta_t$ and $\Delta_{\omega}$ of a signal before discussing any special forms of the uncertainty principle. There is no unique definition of these quantities. With appropriate definitions it can be shown that only the Gaussian signal satisfies the uncertainty principle with equality.

Consider a signal $f(t)$ with Fourier transform $F(\omega)$ satisfying

$$\int_{-\infty}^{\infty}f^2(t)dt=1\quad\textrm{(unit energy)}\\
\int_{-\infty}^{\infty}t|f(t)|^2dt=0\quad\textrm{(centered around }t=0)\\
\int_{-\infty}^{\infty}\omega|F(\omega)|^2d\omega=0\quad\textrm{(centered around }\omega=0)$$

None of these conditions is actually a restriction. They can all be satisfied (for signals with finite energy) by appropriate scaling, translation and modulation.

If we now define time and frequency widths as follows

$$\Delta_t^2=\int_{-\infty}^{\infty}t^2|f(t)|^2dt\\
\Delta_{\omega}^2=\int_{-\infty}^{\infty}\omega^2|F(\omega)|^2d\omega$$

then the uncertainty principle states that

$$\Delta_t^2\Delta_{\omega}^2\ge\frac{\pi}{2}\tag{2.6.2}$$

(if $f(t)$ vanishes faster than $1/\sqrt{t}$ for $t\rightarrow\pm\infty$)

where the inequality is satisfied with equality for the Gaussian signal

$$f(t)=\sqrt{\frac{\alpha}{\pi}}e^{-\alpha t^2}\tag{2.6.3}$$

The equation numbers above correspond to the proof below which is from Wavelets and Subband Coding by Vetterli and Kovacevic (p.80):


Method for finding the best non-increasing approximation function of a given signal

You will find some pointers with the correct wording. There is a lot of statistical literature on monotone or monotonic regression (sometimes called isotonic regression). A more generic term is "shape-constrained estimation".

For instance, a few references:

Constrained statistical inference: inequality, order, and shape restrictions, 2005, Mervyn J. Silvapulle and Pranab Kumar Sen
Polynomial algorithms for isotonic regression, 1997, Victor Chepoi and Daniel Cogneau and Bernard Fichet
Fitting monotonic polynomials to data, 1994, Douglas M. Hawkins

The most common version is a least square fit with piecewise constant or linear functions, see for instance Scikit-learn: Isotonic Regression

  An illustration of the isotonic regression on generated data. The
  isotonic regression finds a non-decreasing approximation of a function
  while minimizing the mean squared error on the training data. The
  benefit of such a model is that it does not assume any form for the
  target function such as linearity. For comparison a linear regression
  is also presented.

With more recent progress in optimization, there are works on fitting under other norms ($\ell_1$, $\ell_\infty$) and contraints on the lower or upper bound of the derivative, which should be either positive or negative, to impose non-increasing or non-decreasing approximations.

[EDIT:2018/12/24] Now, on the question whether the signal should be low-pass or not: if you think that a paramteric function is  too restrictive, a low-pass filter would probably too restrictive as well. Indeed, how do you know that the characteristics or the filter is sufficient? wgy should it be linear? does it cope with the nature of the signal and noise well enough?


I solved this by using a different FFT implementation, one that is using a mixed radix algorithm under the hood ( the project is avalaible here https://github.com/itdaniher/kissfft ). The KissFFT supports both a complex as well as a real FFT so I'll post code for both. The real version is up to 45% faster than the complex one. Note that in order to get the original samples, the results of the IFFT have to be scaled by fft size (fft size doesn't have to be a power of 2 but for the real FFT it must be an even number).

The complex FFT/IFFT

    #include "kiss_fft.h"
    #include "kiss_fftr.h"

        int main(int argc, char *argv[])
            {

                unsigned fs;

                QCoreApplication a(argc, argv);
                cout << "opening ...";
                char *pData= OpenRawA("C:\\test\\250_350_c.raw",fs);

                const unsigned numSamples = 4096;//this doesn't have to be a power of 2
                float* Rx = (float*)pData;

                       kiss_fft_cpx *fIn = new kiss_fft_cpx[numSamples];
                       kiss_fft_cpx *fOut = new kiss_fft_cpx[numSamples] ;
                       kiss_fft_cpx *fIfft = new kiss_fft_cpx[numSamples];
                       //populate the (complex) input buffer
                        for(int i=0; i<numSamples;i++){
                            fIn[i].r = Rx[i];
                            fIn[i].i = 0;
                        }

                        cout << "original data { [20,30,40] } " << fIn[20].r << " " << fIn[30].r << " " << fIn[40].r <<  fIn[50].r << " " << fIn[60].r <<endl;
                        kiss_fft_cfg cfg = kiss_fft_alloc( numSamples ,0 ,0,0 );//configure the FFT function to use a forward FFT

                        kiss_fft( cfg , fIn , fOut );//calculate FFT

                        cfg = kiss_fft_alloc( numSamples ,1 ,0,0 );//configure the FFT function to use an inverse FFT

                        kiss_fft( cfg , fOut , fIfft );//calculate inverse FFT

                       free(cfg);

                cout << "IFFT data { [20,30,40] } " << fIfft[20].r/numSamples << " " << fIfft[30].r/numSamples << " " << fIfft[40].r/numSamples << " " << fIfft[50].r/numSamples << " " << fIfft[60].r/numSamples;
//delete allocated resources
                return a.exec();
            }

The real FFT/IFFT

 int main(int argc, char *argv[])
    {

        unsigned fs;

        QCoreApplication a(argc, argv);
        cout << "opening ...";
         char *pData= OpenRawA("C:\\test\\250_350_c.raw",fs);

        const unsigned numSamples = 4096;
        float* Rx = (float*)pData;

               kiss_fft_scalar *fIn = new kiss_fft_scalar[numSamples];
               kiss_fft_cpx *fOut = new kiss_fft_cpx[numSamples] ;
               kiss_fft_scalar *fIfft = new kiss_fft_scalar[numSamples];
               //populate the input buffer
                for(int i=0; i<numSamples;i++){
                    fIn[i] = Rx[i];

                }
                cout << "original data { [20,30,40] } " << fIn[20] << " " << fIn[30] << " " << fIn[40] <<  fIn[50] << " " << fIn[60] <<endl;
                kiss_fftr_cfg cfg = kiss_fftr_alloc( numSamples ,0 ,0,0 );//configure the FFT function to use a forward FFT

                kiss_fftr( cfg , fIn , fOut );//calculate FFT

                cfg = kiss_fftr_alloc( numSamples ,1 ,0,0 );//configure the FFT function to use an inverse FFT

                kiss_fftri( cfg , fOut , fIfft );//calculate inverse FFT

               free(cfg);

        cout << "IFFT data { [20,30,40] } " << fIfft[20]/numSamples << " " << fIfft[30]/numSamples << " " << fIfft[40]/numSamples << " " << fIfft[50]/numSamples << " " << fIfft[60]/numSamples;
//delete allocated resources
        return a.exec();
    }


Any spectrum that is between the two DFT bin centers ends up getting represented, in various proportions, in all the DFT bins (in the rectangular window case, the proportions decay as in a Sinc function, roughly by 1/(i-k)).  Those proportions might be needed for your synthesis to be an accurate enough reconstruction.

Using 2 bins alone usually creates a beating effect at a low beat frequency.


Receiver function, frequency domain deconvolution not giving logic results

The "spikes" that you see are coming from a "half-done" inverse Discrete Fourier Transform (DFT).
The main thing to note is that when you specify a filter in the frequency domain (as it is hapenning here), that filter has to "fit" the DFT spectrum ordering too.
You conclude (correctly) "I'm splitting the transforms range in the middle since the real input will produce symmetric outputs. " but that is not what is expressed in the code.
If you use numpy.fft.fft() then you have to specify a symmetric filter. Which is clearly not what gauss will be, because it is defined over freq whose linspace goes from 1 (why one?) to numpoints//2. If you plot this, you will see that it is only half the Gaussian and it has to be extended with its mirror image (since you are using numpy.fft.fft()) to include the negative frequencies part.
What you describe (that you only specify half the spectrum because the rest is symmetric) is what the pair numpy.fft.rfft(), numpy.fft.irfft() implement. But if you choose this way of doing it, you have to get rid of numpoints//2 because it is already taken care for you.
Having done this, we then turn to the "cut off" (or equivalent cut off) of the filter. The formulas in Langston79 are in angular frequency and the "standard deviation" (the width of the Gaussian filter) is specified at some $4 \alpha$, which is something like 6Hz (for an $a=1.5$). This has to be defined appropriately on the discrete frequency axis. If you have (for example) an $Fs=1000Hz$ and you are looking for 6 Hz in a 500 point signal, then after the DFT, each bin is $\frac{Fs}{N_{FFT}}$ Hz wide (or 2 Hz in this case). So your discrete $\alpha$ (let's call it $\alpha_d$) is now $3$. If your signal is very long (which I think it is, here, just looking at the duration of those signals), the $\alpha_d$ will work out much higher of course, because the resolution of the DFT is now higher (to still achieve the desired $4 \alpha$).
Hope this helps.


DFT (FFT) of a Real Even Function Doesn't Yield Real Only DFT Signal

When you work with DFT you need to remember one big assumption of the DFT - The samples are part of periodic signal.

For example, assume your signal in on the grid vX = [-3:3].
Let's say it is the simplest symmetric function: vY = vX .^ 2.
If you apply fft on it the result won't be real.

Why is that?
Because the DFT assumes the input is periodic.
Hence the input signal is something like [vY, vY, vY, ..., vY].

Moreover, for the DFT the indexing is [0, 1, ..., N - 1] which in this case is [0:6]. Now if we do periodic continuation, what will you get on the grid of [-7:6]?

figure();
plot([-7:6], [vY, vY]);

It won't be a be symmetric around 0 (Of course look at [-6, 6]).

Let's define vYY = vY(1:6).
If you apply fft(vYY) the result is pure real.
Why is that?
Because for the DFT it sits on the grid [0:5].
If you plot a periodic continuation of it:

figure();
plot([-6:5], [vYY, vYY]);

Look on [-5, 5], now it is symmetric.
Hence its DFT is pure real.

MATLAB Code

% The Grid
vX = [-3:3];
vY = vX .^ 2; %<! Even Function on the Grid

figure();
plot(vX, vY);
title('Figure 001');

vYDft = fft(vY);
norm(imag(vYDft), 'inf')

norm(imag(fft(vY(1:6))), 'inf')

figure();
plot([-7:6], [vY, vY])
title('Figure 002');

figure();
plot([-6:5], [vY(1:6), vY(1:6)]);
title('Figure 003');

Output:

  ans = 7.2678
  ans = 0

Figure 001

Figure 002

Figure 003


Addition of two signals in polar coordinates

Mixing is not the correct terminology for what you want to do.  Mixing is multiplication and you want additive noise.  Convert back to I/Q, add complex noise, then convert back to mag/phase!


Measuring the "frequency stability" of a signal

Your technique might work, provided you window and so on correctly, but it does not really strike me as the right approach for pitch vibrato for a lot of reasons, mostly that the change in pitch is very small compared to other things which could change the "width" of your spectrum (eg, timbre).

I think you'd be better off actually tracking the fundamental pitch in the time domain using a technique like YIN and seeing how much it varies over time.


I believe the short answer is no.

The paper is behind a paywall so I didn't look at it.  In the abstract is the following blurb: " use interpolation techniques to resample the signal at constant angular increments, followed by a common discrete Fourier transform to shift from the angular domain to the order domain."  In which they are describing other techniques.

Looking at your equations, that is precisely what needs to be done to $\theta[n]$ in order to convert it into a regular DFT, suitable for a FFT.

The FFT exploits symmetries in the even spacing of the exponents in order to save computations.  With a varying $\theta[n]$, those symmetries aren't there so you aren't going to be able to implement a FFT algorithm.

You don't tell us what $x[n]$ is either.

How large is N?

Hope this helps at least a little.

Ced


What is Typical RTP Packet size for H.264 video transmission

RTP Payload Format for H.264 Video

   0                   1                   2                   3
   0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  |V=2|P|X|  CC   |M|     PT      |       sequence number         |
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  |                           timestamp                           |
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  |           synchronization source (SSRC) identifier            |
  +=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
  |            contributing source (CSRC) identifiers             |
  |                             ....                              |
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

    RTP header according to RFC 3550

The RTP header information to be set according to this RTP payload
   format is set as follows:

Marker bit (M): 1 bit
      Set for the very last packet of the access unit indicated by the
      RTP timestamp, in line with the normal use of the M bit in video
      formats, to allow an efficient playout buffer handling.  For
      aggregation packets (STAP and MTAP), the marker bit in the RTP
      header MUST be set to the value that the marker bit of the last
      NAL unit of the aggregation packet would have been if it were
      transported in its own RTP packet.  Decoders MAY use this bit as
      an early indication of the last packet of an access unit, but MUST
      NOT rely on this property.

     Note: Only one M bit is associated with an
     aggregation packet carrying multiple NAL units.  Thus, if a
     gateway has re-packetized an aggregation packet into several
     packets, it cannot reliably set the M bit of those packets.

Payload type (PT): 7 bits
      The assignment of an RTP payload type for this new packet format
      is outside the scope of this document and will not be specified
      here.  The assignment of a payload type has to be performed either
      through the profile used or in a dynamic way.

Sequence number (SN): 16 bits
      Set and used in accordance with RFC 3550.  For the single NALU and
      non-interleaved packetization mode, the sequence number is used to
      determine decoding order for the NALU.

Timestamp: 32 bits
      The RTP timestamp is set to the sampling timestamp of the content.
      A 90 kHz clock rate MUST be used.

H.264 also implements B-Frames. This type of frame fills in information between P-Frames. This means that the B-Frame will need to be held until the next P-frame arrives before the information can be used. B-Frames are not used in all modes of H.264. The encoder decides what type of frame is best suited. There are typically more P-frames than I-frames. Lab analysis has shown TelePresence I-frames to generally be 64K wide (50 packets @ 1316 bytes), while P- frames average 8K wide (9 packets at 900 bytes).

Hope this helps


Extracting narrow-band ZigBee signals(4 MHz) from a wide-band WiFi signal(20 MHz)?

Put band pass filters (5 MHz wide) … Re-sample the chunks to 4 MHz.

Don't do that! If you reduce the sampling rate to 4 MS/s, you need to filter to 4 MHz bandwidth anyway. So you could instead just use 4 MHz wide filters and get rid of the resampling filter.


No. It means 1 out of 5 symbols was not able to be decoded. Beyond that it depends on your encoding scheme and criteria for deciding if that packet is usable. For example what if you used an encoding scheme with BPSK symbols such that 1 1 1 1 1 = "1" and -1 -1 -1 -1 -1 = "0"?   (You could replace those with other codes as well 1 1 -1 -1 1 = "1" and -1 -1 1 1 -1 = "0" for example). Typically the packet is longer and error correction encoding is used so that individual symbol errors would not cause the entire packet to be lost.


Notes or books to quickly get started on the basics of signal processing or wireless communication

One of the most readable (imho) authors on the communication systems is Simon Haykin and the following book from him would probably address most of the issues you would encounter in a wireless communication system analysis.

Modern Wireless Communication Systems

Similar known authors do have related books, but I guess most compact academic begining would be this one.

It's not a practical book by any comparison, and therefore does not discuss any bit of a practical application of any of the mathematical analysis it intruduces...

If this book contains advanced materials then you may consider reading its pre-requisites from the same author:

Introduction to Analog and Digital Communications

Furthermore you can enhance your Matlab applications on the fundamental communication systems by the following book from the author Proakis.

Contemporary Communication Systems using Matlab


Lane detection and tracking using B-snake?

Following my comment I suggest you a "dumb" heuristic you could try : Consider your images row-wise. Extract a few features from each row, these could be average RGB value, RGB covariance matrix, entropy etc ... Use a simple clustering algorithm to separate your image into 2 parts : the upper one and the lower one. Let's assume the heuristic has been successful and the lower part is mostly composed of road. We'll just consider this part and run a clustering algorithm on the RGB values. We keep the number of cluster low and we want to take the cluster which is the closest from grey as the road cluster. Once you have a not too bad road segmentation, play around a bit with edge detection to extract only the white marking.


How long does it take for snd_pcm_writei() (ALSA) to send data?

Please correct me if I'm wrong.

snd_pcm_hw_params_get_period_size() gives the number of times per second a chunk of data is sent, i.e. in Hz.

snd_pcm_hw_params_get_period_time() gives the number of microseconds to send such a chunk.

E.g. for my example, it takes in reality around 10 times per second, and these calls return to me 4410 Hz and 10000 usec, respectively. Hence, this appears correct.

However, it still eludes me how the buffer size plays into this. If anyone who happens to be experienced with ALSA would like to add that, that would be awesome.

Edit:

I believe it works like this:

In the same example, 4410 frames will be required for each period. This makes it 8820 samples. If you're using S16_LE, i.e. short data types, then each frame is 2 byte, hence the buffer should be 17640 bytes.

So, that is the buffer which is treated during each period.

However, there is another important concept in ALSA, the hardware buffer, which can be quite large and introduce a latency when copying it---Hence why ALSA uses the concept of periods at all.

Reference: https://www.linuxjournal.com/article/6735

  " The size of the buffer can be programmed by ALSA library calls. The
  buffer can be quite large, and transferring it in one operation could
  result in unacceptable delays, called latency. To solve this, ALSA
  splits the buffer up into a series of periods (called fragments in
  OSS/Free) and transfers the data in units of a period.

  A period stores frames, each of which contains the samples captured at
  one point in time. For a stereo device, the frame would contain
  samples for two channels. Figure 1 illustrates the breakdown of a
  buffer into periods, frames and samples with some hypothetical values.
  Here, left and right channel information is stored alternately within
  a frame; this is called interleaved mode. A non-interleaved mode,
  where all the sample data for one channel is stored followed by the
  data for the next channel, also is supported. "


There seem to be four questions here. I can't give complete answers for each of them but let me address them in turn.

  What exactly is one looking for in audio data to find different instruments or energy sources?

There's not really a single answer to this. There's lots of research into many closely related but ultimately different problems. In general instruments vary in the way their spectral properties vary over time. Tracking the changes of all the harmonic and inharmonic components of a note amongst other sounds is no small task.

  So the FFT can extract frequencies, but how can e.g. percussion be detected alongside e.g. a bass line?

Generally you'd use a Short-time Fourier transform so that there's time data as well as the spectral information. From here you've got a few problems to think about:

Source separation: look up techniques such as ICA, NMF and CASA.This will divide up the energy between the detected sources allowing you to detect the separate instruments.

Fundamental frequency detection: I recommend summary auto-correlation of frequency bands rather than purely time or frequency domain methods.

also Source identification (which I don't know much about).

  Would the frequency bins be checked into a frequency mapping of instruments? This would lead to collisions in possible instruments.

You need to apportion the frequency information not just assign it all to one instrument.

  Are the frequencies of a FFT output directly mappable to pitch frequency charts or should there be a conversion first?

The fundamental frequency, which is given on a pitch-frequency chart, may not actually present in a note ('the missing fundamental problem'). So a fundamental frequency detection algorithm such as the one I suggest above will be needed before you can use the chart.


A first suggestion (without being sure to understand the question).

After sorting on X, try a line fit in a short sliding window (say 4 units long) and observe the correlation coefficient.

In smooth areas, the fit will be excellent. It will degrade at corners. Looking at the plot of the correlation, you should be able to set a suitable detection threshold.


The term delay line comes from (analog) telecommunications where timing between different signals needed to be adjusted by delaying one of them. For example, the chrominance and luminance channels in a cathode-ray-tube colour television.

The reason it's called a line is because that's all it really is: a wire. The website linked to above has the following diagram of one.

In the discrete (digital) domain, I've more usually heard delay lines called tapped delay lines.  That's because each sample is a delay, and to form an FIR filter, we need access (to tap into) the signal at different delays.  The line of $z^{-1}$ boxes in the diagram below is what I would call a tapped delay line.


A narrow-band filter may remove some of the bandwidth of the signal containing timing information of the signal's envelope modulation.  The narrower the filter, the more any tone burst will be smeared out in time.  Upsampling after the timing information is destroyed won't help much.  Using a wider bandwidth signal, channel, and filter, with perhaps a higher sampling rate to support that bandwidth, will help improve time estimation accuracy.

A chirp, or other large frequency modulation is one common method used to increase a signal's bandwidth to allow better timing estimation.

Also note that any real signal with a frequency near either DC or Fs/2 can face strong interference from it's negative frequency conjugate image in the frequency domain, as well as any of the images sidebands, which can also greatly distort the signal envelope and the envelopes timing information in the received sampled data.  A higher sampling rate can also help solve this problem for high frequency tones.  (And a longer sample window length for near DC signals.)


Avoiding clicks with changing biquad coefficients

Another possibility is to send the coefficients as poles and zeros and interpolate/smooth in the z-plane instead of the coefficient domain. This is always guaranteed to be stable. You can either interpolate in real/imaginary presentation or, if you want to get fancy, in magnitude and phase. The first one is a lot easier, since the phase interpolation has to handle wrap around.

Coefficient calculation from real and imaginary part is straight forward and suitable for real time. Only exception may be if you have a dual real pole/zero instead of a complex conjugated pair.


Restricted Isometry Property (RIP) in Compressive Sensing

The restricted isometry property states that:
$$\begin{equation}
  (1-\delta_S)||x||_2^2 \le ||A x||_2^2 \le (1+\delta_S)||x||_2^2
\end{equation}$$
for any $S$-sparse vector $x$. The restricted isometry constant is $\delta_S$, $0 < \delta_S < 1$.

This means that the matrix $A$ is guaranteed to only change the length of any vector $x$ "very little" as long as the vector $x$ is at least $S$-sparse (has at most $S$ non-zero coefficients).

Suppose we have arbitrary $\frac{S}{2}$-sparse vectors $x$. In order to be able to reconstruct such vectors in general, from measurements taken as $y = A x$, we need to be sure that it is possible to distinguish between measurements $y_1 = A x_1$ and $y_2 = A x_2$ of any two such vectors. If $y_1 = y_2$ for any two such vectors $x_1$ and $x_2$, we would not be able to distinguish them and reconstruct them unambigously. Therefore we need to make sure that measurements of any two $\frac{S}{2}$-sparse vectors are "sufficiently different".

If we calculate the difference between any two $\frac{S}{2}$-sparse vectors, their difference can be at most $S$-sparse. So for reconstructing any $\frac{S}{2}$-sparse vector correctly from measurements taken with $A$, the restricted isometry property quantifies how well $A$ lets us do that (the smaller $\delta_S$, the better).

For an early introduction to compressed sensing and the restricted isometry property (and other concepts), see Candès & Wakin, 2008.


Has anyone used Siglib? Answer a question on the FFTr function?

This is the frequency normalized by the sampling frequency, i.e. in your case the frequency in Hz is

$$f=0.15625\cdot f_s$$

where $f_s$ is the sampling frequency in Hz. This can also be found in the SigLib reference manual (p. 251).


Why oversample data before using raised cosine filter?

In digital communications, you can transmit a sequence of numbers using a train of pulses that don't interfere with each other. A given raised cosine pulse $p(t)$ has zero crossings every $T$ seconds around its peak (for a certain $T$). Let's say you want to transmit the number sequence $\lbrace 3, -5 \rbrace$. You can do this with the 2-pulse train $$s(t)=3p(t)-5p(t-T).$$ If you sample $s(t)$ at $t=0$ and $t=T$, you'll get the values 3 and -5.

Your question is related to implementing this idea in Matlab. To do this, you have to solve two different problems: how to create the pulse $p(t)$, and then how to create the train pulse $s(t)$.

Let's start with the pulse $p(t)$. Since you specify Fd=1, then the number of samples in the pulse between the peak and its zero-crossings is given by over. In your case, you can see that the peak occurs at sample 49 and the first zero crossing at sample 65, or 16 samples apart, which is the value of variable over. So, you have defined a pulse $p(t)$ with $T$ equal to 16 sample intervals.

In order to generate the pulse train $s(t)$, you start with the sequence sym. Then you insert 15 zeros between the elements of sym, and finally calculate the convolution of the spaced sequence and the pulse. The result of the convolution will be the addition of pulses $p(t-T)$, scaled with the values of sym.

Let us return to the example of transmitting the sequence $\lbrace 3, -5 \rbrace$. You start with

sym = [3 -5];

from which you generate sig:

sig = [3 0 ... 0 -5 0 ... 0];

The convolution of sig and pulse is the addition of 3*pulse starting at sample 1, and -5*pulse starting at sample 16. The first pulse has a peak equal to 3 at sample number 49, and the second pulse has a peak equal to -5 at sample number 65. The pulses interfere with each other at most sample times but, crucially, the second pulse is zero at sample 50, and the first is zero at sample 66. So, at those precise times, they don't interfere with each other. By sampling the convolution at those times, the receiver can recover the transmitted sequence.

Try this simple code in Matlab:

p=rcosine(1,16,'normal',0.35);
s=[3,-5];
sig=kron(s,[1 zeros(1,15)]);
ss=filter(sig,1,x);
plot(ss);

Using the data cursor, you can verify what I said above. You can also print those values with:

[ss(49) ss(65)]


Any good C# library for signal processing?

I can recommend these libraries:

BASS - also has lot's of functions. Although it is not written for C#, they provide you with a wrapper (Bass.Net). Keep in mind that if you want to use this library in commercial project then you have to pay for the license.

NAudio - is a mature project with lots of functions and support for variety of codecs.

cscore - rather new project with some nice capabilities, especially for visualizations, and some codecs that NAudio is missing. You might want to give it a try at some point, or even contribute to it.


There are soo many projects, field of applications, platforms, architectures and tools to mention about but let me put here a few.

1- Begin with Audio projects, as probably every DSPeng. did so.
There's an open source program AudaCity which involves very rich set of audio functions implemented. So write your own, add your own, develop your own.

Again, audio dsp effects is an area rich for new and better algorithms. Try developing your own.

Digital audio synthesizers is another rich field with immense set of applications.

2- Do image video processing applications. If you are into data compression as well I highly suggest you FFMPEG project. Get into such a project and develop your own applications. There are many filters, effects, algorithms etc that are implemented and you can definietely add your own.

All these require nothing but a PC to develop and implement.

3- In addition, you can search hardware based projects too but they will require a more elaborate set-up. Try implementing some embedded systems using development boards. ADI, ST, TI, NI, NXP provide very convenient, versatile application development boards using their own microcontrollers, microprocessors and DSP procesors, multimedia processors and sensors.

4- Look for some SDR (software defined radio) for most affordable digital communications applications. RTL-SDR is its home. Get a usb-sdr dongle (a hack of DVB-T receiver) and write your own pure C-code decoding algorithm there.

5- A raspberry pi 3 with a pi-camera can provide you a real opportunity to deal with raw camera processing. You can develop a lot of video coding algorithms there.


So, it's generally right that a receiver adds noise due to physics (Johnson-Nyquist noise). So, yes, in any system, at least part of the noise in the received signal is caused by the receiver.

There's also noise that is background noise from e.g. cosmic sources that happen neither in the receiver nor the transmitter (and also aren't other transmitters interfering), but those sources are typically far less important than receiver noise – at least in the microwave bands, where modern comms happen; in LF/MF/HF bands, background noise is very observable.


You won't know the exact channel, instead you'd use the channel estimates.
The basic idea behind Non Orthogonal Multiple Access (NOMA) is that you are receiving multiple overlapping signals, decoding the strongest signal, subtracting it, decode the next strongest signal, and so on. With that being said, NOMA and the associated successive interference cancellation (SIC) receiver requires a high enough SNR so as to be able to reliably decode all of the signals. For this high SNR case, it is reasonable to assume that the channel estimates you have will be good enough for the processing.


Cross Correlation, how can any signals except the trivial cases be uncorrelated?

Tl;DR version:  You are not missing anything; finite-duration signals cannot be uncorrelated signals.

If the crosscorrelation function $x\star y$ of $x$ and $y$ is zero everywhere, then the Fourier transform of $x\star y$, which is $X(f)Y^*(f)$ (or $X^*(f)Y(f)$ for left-handed folks), must also have value $0$ for all $f$. But, finite-duration signals have Fourier transforms with support $(-\infty, \infty)$
and so the product $X(f)Y^*(f)$ cannot be zero everywhere.

It is of course possible to have uncorrelated signals, but these have Fourier transforms with finite (and non-overlapping) support, and so (by duality) must have infinite duration.


State space equation from differential equation

Assuming that the inputs are $f_1$ and $f_2$,

$$\dot q = 0 \cdot q + \begin{bmatrix} \quad\left(\frac{T}{T+1}\right)\\ - \left(\frac{1}{T+1}\right)\end{bmatrix}^\top \begin{bmatrix} f_1\\ f_2\end{bmatrix}$$

Since you have not specified the outputs, the output equation cannot be written.


BIBO stability of $y(t)=\int_{-\infty}^{t}{x(\tau)d\tau}$

Since this is most likely homework, here is a hint.

Write the integral you have displayed in the form $\int_{-\infty}^\infty x(\tau)h(t-\tau) d\tau$ where you get to choose what the function $h(\cdot)$ is to make it all work out. Then, $h(t)$ is the impulse response of the LTI system, Do you know the criterion for BIBO stability of an LTI system in terms of $h(t)$? Does your $h(t)$ satisfy the criterion?

Edit in response to OP's comment and definition of BIBO stability:

Your understanding of the definition of BIBO stability is faulty.

  A system with input $x(t)$ and output $y(t)$ is said to be bounded-input bounded-output (BIBO) stable if it has the property that whenever $x(t)$ is bounded (that is, there exists a (finite) real number $M$ such that $|x(t)| \leq M$ for all real numbers $t$), then the output is also
  bounded (that is, there exists a (finite) real number $N$ (possibly dependent on $M$) such that $|y(t)| \leq N$ for all real numbers $t$).

Note that while $N$ might depend on $M$, it does not depend on $t$, the same $N$ must be an upper bound on $|y(t)|$ for all $t$.

In your case, with $y(t) = \int_{-\infty}^t x(\tau) d\tau$, what happens if $x(t) = u(t)$ for all real numbers $t$ (step input)? This is a bounded input. Does there exist a finite number $N$ such that $|y(t)| \leq N$ for all real numbers $t$?


Signal processing techniques for an accelerometer signal?

I'm wondering why the STFT pops out. To me, wouldn't a simple threshold on the signal itself or on its envelope do better / just as well, after removal of the g offset?
Once you decide what "measure" is best to detect your event, you can apply the work of Basseville and Nikiforov, that I answered here.

The classic reference for that problem is Detection of Abrupt Changes - Theory and Application by Basseville and Nikiforov.  The whole book is available as a PDF download.
My recommendation is that you read Chapter 2.2 on the CUSUM (cumulative sum) algorithm.


How to calculate step response for $y''(t) - y(t) = x'(t) - x(t)$ in time domain?

It is indeed possible to convert the system to an equivalent first order system, but let's first directly solve the given problem. The system's step response satisfies

$$y''(t)-y(t)=\delta(t)-u(t)\tag{1}$$

This implies that we can compute the step response of the given system by considering another system

$$y''(t)-y(t)=x(t)\tag{2}$$

and compute its impulse response ($x(t)=\delta(t)$) and its step response ($x(t)=u(t)$), and obtain the step response $(1)$ as the difference of the impulse response and the step response of the new system $(2)$.

For the step response we need to solve

$$y''(t)-y(t)=u(t)\tag{3}$$

Assuming the system is initially at rest, we get initial conditions $y(0)=y'(0)=0$. The homogeneous solution of $(3)$ is $y_h(t)=Ae^{t}+Be^{-t}$, $t>0$, whereas a particular solution is $y_p(t)=-1$, $t>0$. The complete solution is then

$$y(t)=y_h(t)+y_p(t)=Ae^{t}+Be^{-t}-1,\quad t>0\tag{4}$$

The constants $A$ and $B$ are determined by the initial conditions:

$$\left.\begin{align}y(0)&=A+B-1=0\\y'(0)&=A-B=0\end{align}\right\}\Rightarrow A=B=\frac12\tag{5}$$

Consequently, the step response of system $(2)$ is

$$y(t)=\frac12\left(e^t+e^{-t}\right)-1,\quad t>0\tag{6}$$

Its impulse response satisfies

$$y''(t)-y(t)=\delta(t)\tag{7}$$

and it can be obtained by differentiating the step response $(6)$:

$$y(t)=\frac12\left(e^t-e^{-t}\right),\quad t>0\tag{8}$$

Finally, according to $(1)$ and $(2)$, the step response of the original system is the difference between the impulse response and the step response given in $(8)$ and $(6)$:

$$y(t)=\frac12\left(e^t-e^{-t}\right)-\frac12\left(e^t+e^{-t}\right)+1=1-e^{-t},\quad t>0\tag{9}$$

Another approach would be to see that the original system can also be represented by a first-order system. Note that

$$y''(t)-y(t)=(y(t)+y'(t))'-(y(t)+y'(t))\tag{10}$$

So if $y''(t)-y(t)=x'(t)-x(t)$, it follows from $(10)$ that

$$y(t)+y'(t)=x(t)\tag{11}$$

With $x(t)=u(t)$ you obtain from $(11)$ the same step response as in $(9)$:

$$\begin{align}y_h(t)&=Ce^{-t}, &t>0\\
y_p(t)&=1, &t>0
\end{align}\\y(t)=y_h(t)+y_p(t)=Ce^{-t}+1,\quad t>0$$

From $y(0)=0$ we get $C=-1$, and the step response becomes

$$y(t)=1-e^{-t},\quad t>0$$

just like in $(9)$.

As a final check it is useful to see what happens when you solve the problem using the Laplace transform. The original difference equation transforms to

$$Y(s)(s^2-1)=X(s)(s-1)$$

which gives for the transfer function

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

This shows that you have a pole-zero cancellation, and the system is actually a first-order system, just as described by Eq. $(11)$. For $x(t)=u(t)$ you get $X(s)=1/s$, and, consequently, for the Laplace transform of the step response

$$Y(s)=\frac{1}{s}\frac{1}{s+1}=\frac{1}{s}-\frac{1}{s+1}$$

which corresponds to the time domain function

$$y(t)=u(t)-e^{-t}u(t)=(1-e^{-t})u(t)$$

which is the same as $(9)$.


DC, by definition, is the non-varying part of the signal. To find the DC component you simply average the entire signal.

  How can I extract the varying DC wave from the signal?

This sounds more like you want to extract the low frequency part of the signal. The most basic filter for this is a moving average filter. However, I would recommend a smoother filter, such as a Gaussian filter. Start with $\sigma$ about 1/6 the period of the low frequency oscillations you are interested in and modify to suit.


How to accurately model wideband RF reflections?

So, to give you something to read up on first, the channel you describe is a Rician or Rayleigh channel, depending on whether you have a dominant line-of-sight path or not.

So, as a first approach, to delay something in time, you don't have to shift it by a whole sample – you can also do it in frequency domain, by DFT'ing your signal, multiplying it with a $e^{j\frac{\Delta t}{f_s} f}$, and then IDFT it back. That's how I've built all my radar simulations.

  I'm assuming that each reflector can be modeled by the distance to the transmitter and the receiver, and some real-valued attenuation constant.

A radar target typically also phase-shifts the electromagnetic wave – on perpendicularly hitting a metal surface, by 180°, but differently for more complex surfaces. In reality, we often model that stochastically. What is appropriate for you will very much depend on what you simulate. A model for scattering of reflections e.g. on leaves will not be very useful for modeling reflections e.g. inside an oil tank, and vice versa. Both are use cases for radar.

For many scenarios, there's "standard" models, such as the COST-Hata urban model, the ITU indoor models and many, many more. It really depends on which environment you're considering, and also, which group of researchers you ask.

Many of these models are defined for what I call the "easy RF" ranges below ~5.8 GHz. Since radar, for Doppler and wavelength reasons, often operates above, things often aren't too well represented by these models (i.e. they simply don't apply to common radar frequencies > 8GHz). This has two interesting aspects:

The higher your frequency, the more quasi-optical your models get for most useful spatial scales, and
Since millimeter waves are very en vogue right now (with all the 5G research going on), you might find very interesting RF propagation models when looking at very recent papers. As an example for a very practical, hands-on-testing based propagation/reflection description, try Simić, Perpinias and Petrova: 60 GHz Outdoor Urban Measurement Study
of the Feasibility of Multi-Gbps
mm-Wave Cellular Networks. Downside of these considerations is always that, albeit multipath does play an important role for these researchers, the individual phase really doesn't matter (only the interference with other paths does).


FFTW factorizes by testing against a table of primes.  This is fast due to the fact that the size of the table is limited because the size of computer memory in which the FFT data is to be computed is itself limited to a very tiny number in mathematical terms (less than a few petabits). You only need up to the square root of that number.


FIR filter Gain of 0dB at Passband Frequency

FIR coefficients are $h[n]$, the same as the impulse response.  there are $N$ non-zero taps.

$$ h[n] = 0 \qquad \text{for } n<0, n \ge N $$

frequency response is

$$\begin{align}
H(e^{j \omega}) &= \sum\limits_{n=-\infty}^{+\infty} h[n] e^{-j \omega n} \\
                &= \sum\limits_{n=0}^{N-1} h[n] e^{-j \omega n} \\
                &= \sum\limits_{n=0}^{N-1} h[n] \big( \cos(\omega n) - j \sin(\omega n) \big) \\
\end{align}$$

the magnitude-squared is

$$\begin{align}
\Big| H(e^{j \omega}) \Big|^2 &= \left| \sum\limits_{n=0}^{N-1} h[n] e^{-j \omega n} \right|^2 \\
                &= \left| \sum\limits_{n=0}^{N-1} h[n] \big( \cos(\omega n) - j \sin(\omega n) \big) \right|^2 \\
                &= \left| \sum\limits_{n=0}^{N-1} h[n] \cos(\omega n) - j \sum\limits_{n=0}^{N-1} h[n] \sin(\omega n) \right|^2 \\
                &= \left( \sum\limits_{n=0}^{N-1} h[n] \cos(\omega n)  \right)^2 + \left( \sum\limits_{n=0}^{N-1} h[n] \sin(\omega n) \right)^2 \\
\end{align}$$

you must square root that expression to get magnitude of gain

$$ \Big| H(e^{j \omega}) \Big| = \sqrt{\left( \sum\limits_{n=0}^{N-1} h[n] \cos(\omega n)  \right)^2 + \left( \sum\limits_{n=0}^{N-1} h[n] \sin(\omega n) \right)^2 }  $$

evaluate at at your known PassbandFreq $\omega_0$ (normalized angular frequency) and divide each $h[n]$ by that gain magnitude:

$$ h[n] \leftarrow \frac{h[n]}{\ \big| H(e^{j \omega_0}) \big| \ } $$


i dunno what a "Lipshitz filter" is but i do know who Stanley Lipshitz is.  now, i don't know for sure what Stanley would say to you, but the best advice i got from him is that it does little good to color the additive dither if you are also using error feedback.  if there is no noise shaping feedback, probably the simplest way to get triangular PDF dither is to generate good, independent, rectangular PDF random numbers $v[n]$ and run that through a differentiator.

$$ d[n] = v[n] - v[n-1] $$

the result is high-pass in the spectrum and triangular PDF.  if $\Delta$ is the step size of the quantization, $\frac{\Delta^2}{12}$ is the variance (or power) of raw quantizer, but you have two more of them due to your dither and $\frac{\Delta^2}{4}$ is the total power due to the entire dithered quantizer.  a 4.77 dB increase in quantization error, but the quantization is fully decoupled from both the mean and variance of the input signal.  although this is the same quantization power that you would get from adding two independent random numbers, i think you would prefer that much of the power is steered to higher frequencies (but not all of it, the spectrum of the raw quantizer is still flat).

if you're planning on error feedback, i think the best advice is either use white, triangular PDF dither (requires adding two independently generated rectangular PDF random numbers) and error feedback (designing a good feedback filter takes a bit).

or no dither (if the quantization step size is very small) and error feedback.  the simplest no-trivial example of this is "fraction saving" where you simply hack off the low-order bits (the floor() function) and add those bits back in (zero-extended on the left) to the next sample before that sample is quantized (with the floor() function).  it has infinite S/N ratio at DC and totally kills any limit cycle in feedback filters (like biquads) that get the filter stuck on a DC value even when the input goes to zero.  i've done that a lot.


How do we actually obtain samples to form discrete time samples?

Assume you are speaking, those vibrations (sound) are converted into electrical signals(continuous time signal) by Mike. After that, a Sample & hold circuit is placed, which stores the converted electrical signal inside a capacitor in the form of charge (now u get discrete time signal).

Every signal is  sampled at some frequency called sampling frequency. Let the sampling frequency in this case be 8KHz. It means, at time t=0, you sample (read) the electrical signal amplitude (by using sample & hold circuit) and hold it(store it) in a capacitor. You wait till some time and when time t=0.00125 seconds(1/8KHz, remember our sampling frequency is 8KHz), you again sample the electrical signal and store it . You repeat this process every 0.00125 seconds. Now you have got discrete version of your continuous time electrical signal (in fact speech signal).

In the time between every samples(0.00125 seconds), the charge stored in the capacitor is converted to digital binary values using special circuits called ADCs. These converted digital values and discrete signals in binary form(1s and 0s which may be stored or for further processing in a computer).

If you are interested, read about how sample and hold circuit works. You need some electronics background for that.


µ-law is defined for real-valued PCM signals only. You will either need to produce a real-valued signal containing OFDM as passband signal, either

by mixing it to a carrier frequency, or
by modifying your OFDM system to only produce symmetrical spectrum
baseband, or
by u-law coding I and Q separately (but that makes little sense, since a high PAPR happens when both I and Q have a high value at the same time, i.e. coding them separately isn't logical),

but all three options wouldn't be "companding OFDM" anymore, that would be "companding a system which, among other things, involves OFDM".

General remark: Don't forget that you've got channel coding. It's often better to intentionally completely "delete" isolated subcarriers from extreme-peak OFDM symbols than to  "make whole symbols slightly worse", as the former leads to a certain, but isolated, error, whereas the other has a  chance of leading to multiple errors in the same symbol – which would, assuming there's no large interleaver at work, look like a burst error.

I think you'd want to read

X. Wang, T. Tjhung and C. Ng, "Reduction of peak-to-average power ratio of OFDM system using a companding technique", IEEEE Transactions of Broadcasting, vol. 45m, no. 3, Sept. 1999.

on what they did and what their results are.

While you do that, look at Figure 5:

Notice how the second (⛛) curve is the uncompanded one. This looks like for anything but very high SNR, the companded curves are better.

Now, read

A. Mattsson, G. Mendenhall and T. Dittmer, "Comments on reduction of peak-to-average power ratio of OFDM system using a companding technique", IEEE Transaction on Broadcasting, vol. 45, no. 4, December 1999.

and especially this paragraph:

  Which shows that there is a significant increase in power
  for the companded signal. I.e. the SER for the companded
  signal should in our opinion be compared to the SER for
  the non-companded signal amplified by 10.9 dB [for the µ=16 case], 17 dB [for the µ=255 case] etc. It appears that the SNR refered to in Fig. 5 on page 307,
  refers to the signal power of the signal before companding
  and not the actual SNR of the transmitted signal.

In other words: shift the red (uncompanded) curve 11 dB to the left to compare against the µ=16 case, and by 17 dB to the left to compare against the (standard) µ=255 case.

Look how the uncompanded signal (red line) is always better at the same power than the µ=16 companded signal, and how it's better for SNRs > 12 dB for the µ=255 case. I highlighted the amount that uncompanded is better than companded in green for µ=16, and in orange for µ=255.

If your SNR is below 12 dB, you get a SER > 10⁻³, which is a typical design goal for wireless communication; you might want to use a lower modulation than 16-QAM then, which will then typically be an 8-PSK or QPSK. Which is indifferent about  amplitude variation – so instead of companding nonlinearly, you'd just look at the peak in your OFDM symbol and simply either drastically scale a few isolated subcarriers, or the whole symbol slightly. (which you do would depend on how steep the BER curve of your channel decoder is at that point.)

TL;DR: Literature research shows companding only works for low SNR. There, you might not be in a sweet spot for system design, so my expectations for companding to work are low.  Instead, I propose you generally rather design your FEC so that it can live with intentionally broken isolated symbols that you insert on purpose to reduce PAPR.


Relationship between matrix rank and beamforming

Typically, when doing any sort of adaptive bamforming, one needs to invert a (square) (covariance) matrix and it needs to be full rank in order to do that. Actually there are work arounds if it isn't full rank and it doesn't always require a literal inversion, like using rank one updates of QR or Cholesky decomposition.  So, that leads into what is a rank one matrix? If $\mathbf{x}$ is a column matrix, the outer product , $$ \mathbf{X} = \mathbf{x} \mathbf{x}^H
$$ is a rank one matrix. A matrix of rank $N$ is, assuming each $\mathbf{x}$ span a linear space, is a sum of $N$ rank one matrices.


Median filtering is non-linear, and pretty awesome about removing outliers.  You just need to adjust the length of the filter based on the estimate of the frequency of the errored samples.


fir filter - mean delay and snr from difference equation?

I'll give you some hints to help you solve the problem by yourself. Looking at the input-output relation, you can see that the output signal is computed only from a combination of delayed input signal values, and not from delayed output signal values. So your impulse response has a finite length and the filter can be implemented by a non-recursive structure (FIR filter). Using the impulse response, you can write the input-output relation as a finite convolution sum

$$y[n]=\sum_{k=0}^{N-1}x[n-k]h[k]\tag{1}$$

Comparing (1) to your difference equation, you can easily see what the system's impulse response $h[n]$ is. Once you have the impulse response, you can compute the frequency response using an FFT. Plot its magnitude and phase to see how the filter behaves in the frequency domain.

The (group) delay is the negative derivative of the filter's phase response. You will find that the phase is a linear function, so the delay is a constant. You just need to determine the slope of the phase (either analytically or by inspection of the phase plot).

Without any further knowledge of your desired signal and the noise characteristics, you cannot say much about a possible SNR improvement.


Signal response in Frequency Hz from angular frequency

have a look to this example: http://www.mathworks.com/matlabcentral/fileexchange/55822-low-pass-filtering

here you have a signal that can be represented/plot in both frequencies scales:

omegas_axis=[0:2*pi/L:2*pi-2*pi/L];%axis input signal
frec_axis=[0:fs/L:fs-fs/L];%axis input signal

the "omegas_axis" is the angluar frequency.

Hope this helps,

BTS.


What does the normalization step of the Haar wavelet transform represent?

As I understand it, the normalization is because the Haar wavelet conserves energy of the signal. In that, when you take signal from one domain to another, you aren't supposed to add energy to it, (although conceivably you might lose energy).

The normalization is just a way to ensure that the energy of your Haar-transformed signal in the Haar-domain has the exact same energy as your signal in the original domain.

Intuitively speaking, Haar, Fourier, etc, are all just basis transformations, which intuitively just mean that you are looking at the signal in a different way, (technically, through a different set of bases). Therefore, if all you are doing is looking at a signal differently, its energy cannot/should not change.


Separate static echos from acoustic signal (ICA?)

This probably should be a comment, but I didn't earn that option yet. You can try measuring instrument impulse response. If I understand correctly, your system is essentially one-dimensional. Then if you use thin plate or plate of known thickness instead of thick object/material, your output signal will be an approximation of impulse response. Using it as deconvolution kernel will allow you to measure response specific to material, or object.


You're not doing anything conceptually wrong, you're just misunderstanding MATLAB. The code as you provided it doesn't execute properly; here's a fixed script:

% don't need this; MATLAB already knows what j is
%j = sqrt(-1);
N = 10^6; % number of symbols
Es_N0_dB = 12; % ratio of signal energy to noise energy in dB
QP = (2*(rand(1,N)>0.5)-1)+j*(2*(rand(1,N)>0.5)-1); % actual QPSK data
s = (1/sqrt(2))*QP; % normalization of energy to 1
n = 1/sqrt(2)*(randn(1,N)+j*randn(1,N)); % white Gaussian noise, 0 dB variance
y = s + 10^(-Es_N0_dB/20)*n; % recovered noisy data
plot(y, 'b.')
sr=real(s); si=imag(s);
yr=real(y); yi=imag(y);

% it seems like you want to predeclare these since you're indexing them below
errorI = zeros(size(s));
errorQ = zeros(size(s));

for i = 1:N
    if sr(i)*yr(i)<0
        errorI(i)=yr(i);
        if si(i)*yi(i)<0
            errorQ(i)=yi(i);
        end
    end
end

err = errorI+j*errorQ;
% BER is the number of nonzero entries in `err` divided by the number of chances
BER = sum(err ~= 0) / length(y);


Choosing sampling frequency and base interval to eliminate leakage

The key attribute you need to consider is for a DFT not to have leakage a sinusoidal signal must have a whole number of cycles in the sample frame.  Looking at your frequencies, clearly a one second interval satisfies this criteria, but it is not the shortest interval to do so.  The base interval is the shortest one that you are looking for that has a whole number of cycles for each signal.

Notice, this does not determine the sampling frequency.  Your $0th$ sample and your $Nth$ sample have to occur at the same relative point in your repeat pattern.  To get this, once you have found your base interval in continuous time, divide it into $N$ pieces to get your sampling frequency.  Any choice of $N>1$ will work to not have leakage, but if you want the bin indexes to represent the actual frequencies and not aliases you have to select $N$ large enough so that the higher frequency signal is at or below the Nyquist frequency.  Another way to say this is that you have to have at least two samples per cycle.

Followup:

This is the spoonfeeder version so the OP should stop reading here until they've tried to solve it using my clues themselves.

Honestly, I think Fat32's solution is a little overcomplicated, but impressive.

Let's call the two component sinusoidals Tone1 and Tone2.  Tone1 has 1400 cycles in a second and Tone2 has 1450 cycles in a second.  Observe:

$$ \frac{1400}{1450} = \frac{50 \cdot 28}{50 \cdot 29} =  \frac{28}{29} $$

From this we can deduce that it takes Tone1 28 cycles and Tone2 29 cycles to span the base interval of $1/50$ of a second.  Thus any DFT on this interval, no matter the sampling rate, will have non-zero bin values at bin 28 and bin 29 (or their aliased equivalent) and at the corresponding mirror image bins in the negative direction.

If $N \ge 29$ then it will be bins 28 and 29, otherwise you have to take modulo N to figure out the aliased bin index.  In order to be below the Nyquist limit $N \ge 58$.  Note that if either tone falls on the Nyquist bin (or the DC bin), the results may be zero or all the way up to the actual amplitude depending on where the sampling occurs within the sine wave.  For the OP's equation, the result will be zero since the sampling will occur at the zero crossings.  There still won't be any leakage.

For Fat32's first example, a samping rate of 4000Hz was used.  For an interval of $1/50$ of a second it means $ N = 80 $.  As per the first diagram, you can see that bins 28 and 29 are non-zero, and the mirror image bins 52(=80-28) and 51(=80-29) are also non-zero.

Fat32's second example is actually over an interval of 1/25(=15/375) of a second, so it is twice as long as the base interval.  Thus Tone1 has 56 cycles in the frame and Tone2 has 58.  Since N=15 in the example, aliasing does occur.  The bin indexes are then 11(=56 mod 15) and 13(=58 mod 15).  The mirror image bins are then 4(=15-11) and 2(15-13).  Doubling the interval doubles the spacing between the bins.

Side comment:  Fat32's example also shows why I prefer a $1/N$ normalized DFT.  If it were used, the bin magnitudes would be $1/2$ regardless of $N$.

Response to Fat32's comment:

Your approach is still overly cumbersome.  I think it is because you are thinking about the problem backward, and then backsolving.

Step 1. Set Sampling Rate

Step 2. Find N

Step 3. Calculate Interval length (Impicit)

For your first case:

Step 1. $F_s=4000$

Step 2. $N=80$

Step 3. $T_{frame} = 1/50 = .02$ seconds

For your second case:

Step 1. $F_s=375$

Step 2. $N=15$

Step 3. $T_{frame} = 1/25 = .04$ seconds

Clearly your second case frame interval is twice as long as your first case.  Neither your 375 sample rate or the 4001 sample rate will result in an interval that is the .02 second minimum, thus the 28 and 29th bin stipulation is not applicable.

In contrast, I did it in reverse order to yours:

Step 1. Find Interval

Step 2. Set N

Step 3. Calculate Sampling Rate

The minimum frame interval for which both tones have a whole number of cycles is independent of the sampling rate.  Thus, it can be found before any step into the discrete domain needs to be taken.  Finding the minimum value of .02 seconds was done above and I don't think too difficult to understand.  The only intervals that lead to leakage free DFTs are whole integer multiples of the minimum value.  Any non-integer multiplier will lead to non-whole cycle tones.

Once an interval has been chosen, it quite easy to set the number of sample values (N) to any value you want.  It has to be an integer by definition.

The DFT doesn't care about the sampling rate, it only cares about the number of samples.  The sampling rate is calculated from the frame interval duration and the number of samples.

The most important concept is that the bin index (zero based, not MATLAB's stupid one based) corresponds to the number of cycles per frame.  Thus, no matter what N is, as long as N is greater than 29, it will be bin 28 and 29 that have non-zero values for any DFT done on a .02 second interval of this signal.  The value of N will definitely determine where the mirror image bin values are at.  In your second case, since then frame interval is twice as long, the two tones will have frequencies of 56 and 58 cycles per frame respecively.  So those are the bin indexes (or modulus N thereof) that will have non-zero values along with their mirror images.

All possible solutions can then be found by:

Step 1. $T_{frame} = .02 \cdot m $

Step 2. Set an $N$

Step 3. $F_s=N/T_{frame}$

Where $m$ and $N$ are integer values.

Response to Fat32's followup comments:

Comment 1: You are absolutely correct.  Poor wording on my part.  What I meant was "determined" not "constrained".  I have fixed it.

Comment 2: The OP said "The task is to chose the sampling frequency and the base interval of the discrete signal so that there will be no leakage-effect after calculating the DFT."  I think the salient concept behind the question is that you need to have a whole number of cycles of the constituent tones in order not to have leakage in the DFT.  The sampling rate/number of samples is not important in this regard.  I never called your approach wrong, just cumbersome and overly complicated.

Technically, your approach does not work for any $F_s$.  Integer yes, even rational yes, but irrational no.  There are an uncountable infinity worth of values that won't work compared to a mere countable infinity that will.

I think this is poorly worded:  "the resulting N will take more samples than base interval minimum." Perhaps "the resulting N will take a longer interval than the minimum interval."

C3: I did cover the possibility and consequences in my answer of the bin values falling in the DC or Nyquist bins for some choices of N.  Per the criteria of the problem statement, what is important is that there will still not be leakage if this happens.

I am not a MATLAB user, but I am aware of the one based indexing and am also aware that it is the most common criticism of MATLAB.  If it is the tool of choice for you, fine, but arguing what is the best platform quickly devolves into a religious war and is supposed to be off-topic for SE.

SUMMARY:  The problem wasn't stated in either of those ways.  If either additional criteria were given then that would be a discriminant on which solution to choose.  Absent that, I believe my solution demonstrates the underlying concepts clearer and definitely requires fewer calculations.


Contour and area, raw (spatial) and central image moments

Actually, I was surprised how hard was it do deduct a proper definition of contour versus "normal", non-contour moments of an image. After reading a bunch of materials, here come my conclusions.

Firstly, in order to understand moments, and especially the difference and the usage of spatial (what the OP calls "raw"), central, and central normalized moments, I found two very good materials:

(manual) Johannes Kilian: "Simple Image Analysis By Moments"

Excellent manual with simple mathematics. Don't be scared by the integrals - you can read all of them as summations.

Also, it has a small overview on OpenCV functions used to operate with this moments. It's very old material (2001), so the OpenCV manual it is referring to is a bit old, but it still helps.

And than there's the wonderful third chapter, specifying which moment is used to describe which characteristic of a moment.
(image processing blog) Utkarsh: Image Moments

Simple, short and friendly. I found a lot of good material on this blog before.

Disclaimer AI Shack seemed to be offline at some point. Here is the homepage from the AI Shack author, where he talks about this project, so it still seems to be supported. I hope it comes back online soon, but if not maybe it can be tracked through the author's webpage.

Shortly, the spatial moments give information about the object in the image, i.e. related (dependent) on the object position.

The central moments are adjusted for translational invariance, by moving the origin of the "coordinate system" used for calculations to the centroid (center of gravity) of the object in question.

Finally, the central normalized moments are scaled by the area of the object, and are thus scale invariant in addition to translational invariance.

Now for the actual question part: what about contour moments?

The deductions from this part are mostly based on

OpenCV reference manual, v. 2.4.4.0, Chapter 3.5: Structural Analysis and Shape Descriptors
wikipedia page for Green's theorem - following the link from OpenCV manual, the method used when calculating contour moments (I'll explain in a few lines)
Gary Bradski, Adrian Kaehler: "Learning OpenCV: Computer Vision with the OpenCV Library"

(Link is to google books, but the relevant pages are accessible. The section that I linked to, plus the next 2-3 sections are relevant. That's about 3 pages total)

And the most important quotes from those sources:

  The moments of a contour are deﬁned in the same way but computed using the Green’s formula.

(OpenCV reference manual)

  In plane geometry, and in particular, area surveying, Green's theorem can be used to determine the area and centroid of plane figures solely by integrating over the perimeter.

(wiki for Green)

  Moreover, cvContourMoments is now just an alias for cvMoments.

(Bradski Kaehler book)

Based on that, I would deduce that contour moments does not refer to special measures of the object contours, but instead to a particular way to calculate image moments, only using the contour information (instead of pixel information for the whole image).

The difference, in the fundamental case, would be how both are calculated.

My guess would be that the direct implementation would work by pixel-by-pixel summation, directly implementing the formula. The object is expected to be filled.
My guess for the contour moments would be that the image contours are first determined (consult OpenCV manual) and then the Green theorem is applied on the contour data.

That would make the measurements slightly different for real images because the methods would differ in their: sensitivity to: noise, scaling, discretisation (pixel grid instead of continuous image). Also, the speed: calculating using contours is faster than using the direct approach. I would speculate that they would give perfectly equal results for an (idealized) continuous black and white image with no noise.

So, to answer your questions: the moments should be the same (differing because of noise etc). You can use spatial (raw) moments calculated by both methods to determine central moments (that will still describe the same thing).

Further support of this claims is the existence of this article (I only read the abstract, but should be very relevant, and even the abstract is informative) from 1994:

Luren Yang, Fritz Albregtsen: "Fast and Exact Computation of Moments Using Discrete Green's Theorem"

Note about getting the perimeter measure: I think, to get the "perimeter" which is actually just area of the contour, I would calculate the $0^{th}$ moment of the image of the contours of the objects, but treat the contours as a really thin object, instead of as "contours of an object".

All further measurements would of course differ if you used this moment further.


Your mistake is that $\mathcal Z$ does depend on process noise. It's just a bit more obscure than the linear filter. $\mathcal Z$ is the projection of the sigma points through the observation function. The sigma points are a $\mathcal f(P)$. Therefore $\mathcal Z$ is a $\mathbf f(P)$, and therefore related to $\mathbf Q$.


Your code is correct.

This can also be implemented by calling the built-in MATLAB function hamming:

WIN_HAM_v = round(4095*hamming(LUT_length));

You can (should) compare the two methods yourself:

win1 = 4095*(0.54-(0.46*cos(2*pi*n/(LUT_length-1))));
win2 = 4095*hamming(LUT_length);
maxError = max(abs(win1(:)-win2(:)));
% result: 2.2737e-12


Building periodic sequences from finite-support sequences


Reconstruction of a windowed signal (specifically Hamming windowed)

did you ask this at comp.dsp?  appears so.  if you do what hotpaw says, your Hamming window will be

$$ w(t) =
\begin{cases}
0.426 \cdot \cos(\pi t) + 0.5 & \text{if }|t| < 1 \\
0.037 & \text{if }|t| = 1 \\
0 & \text{if }|t| > 1
\end{cases} $$

still a Hamming window, but it overlap-adds to 1 (with 50% overlap), rather than 1.08.


I figured out the problem. But there is one constraint

(1) The scaling factor for the entire matrix has to be the same number to do the process efficiently

So the way correlation is done, is multiplication and accumulation of two matrices, take for example the image matrix below correlated with the signal, and the result is from the valid part of the correlation, as in no zero padding is done

image        signal    scaling constant = 3
1  2  3      9 4
4  5  6      3 4

result[1] = (1*9+2*4+4*3+5*4)/3 = 16.3333
result[2] = (2*9+3*4+5*3+6*4)/3 = 23

If I want to scale before the correlation, I would have to scale the whole correlation window, like so

image              signal    scaling constant = 3
1/3  2/3  3/3      9 4
4/3  5/3  6/3      3 4

result[1] = ( (1/3)*9+(2/3)*4+(4/3)*3+(5/3)*4 ) / 3 = 16.3333
result[2] = ( (2/3)*9+(3/3)*4+(5/3)*3+(6/3)*4 ) / 3 = 23

After scaling like above, I could alternatively go to fourier domain with the image and the signal and do correlation there, do ifft, and I still get the same results.

image              signal    scaling constant = 3
1/3  2/3  3/3      9 4
4/3  5/3  6/3      3 4

image_fourier=fft2(image,4,6);   zero padding
signal_fourier=fft2(signal,4,6);
result=ifft2(image_fourier.*conj(signal_fourier));

NOTE: above is matlab syntax

result[1] = 16.3333
result[2] = 23;

But this will not work if the scaling factor is different for result[1] and result[2].


A phase measurement needs to be referenced to some point.

The standard FFT computes a circular phase relative to the first sample, or the sample at the 2D origin.  For data that is not exactly integer periodic in the FFT aperture (typically true for any rectangularly windowed data), this is at a circular discontinuity in the data vector.  Thus the result is an alternation of FFT phase results to represent this edge discontinuity.

An edge discontinuity can be represented by an odd function, and odd functions show up in the imaginary (or sine) components of a DFT of this discontinuity.

You can move this circular or edge discontinuity away from the FFT’s phase zero reference point by doing an fftshift before performing the FFT.  This moves the phase reference point to the center of your data, which hopefully does not have an unintended discontinuity.  This shift most often results in a more continuous phase spectrum across FFT result bins of windowed data.  However a pre-fftshift produces the same result as post flipping the phase of all the odd numbered FFT result bins without the shift, which may be easier to do.


Why do most digital pll's only use the positive zero crossing?

Getting accurate lock when using both edges or zero-crossings requires that the duty cycles of both the PLL local oscillator and the reference signal match and/or be as stable as the frequency.  If the duty cycles don't match when using both edges, then you would have to filter out some twice per cycle jitter which could well slow down lock speed.


unit step signal as sum of impulses from MIT OCW video

yes there is only one term in the sum when $n \ge 0$
$$ \delta[m]=\begin{cases}1 \quad m=0 \\
0 \quad \text{otherwise} \end{cases}$$

$$ \begin{align}\sum_{m=-\infty}^n \delta[m] &= \delta[-\infty]+ \dots +\delta[0] + \dots + \delta[n] \\
& = 0 + \dots + 1 + \dots + 0  \\
&= 1  \; \text{for} \; n\ge 0
\end{align}$$
The sum holds the value of 1 for $n\ge0$. His figure shows the terms inside the sum, not the sum of terms.

The definition of $u[n]$ is
$$
u[n]=\begin{cases} 1  \quad \text{for}\; n \ge 0 \\
0  \quad \text{for}\; n<0 \end{cases} $$
They are equivalent.  Defining $u[n]$ function as a sum of delta functions may seem strange but I'm sure he has a "trick" to show that will be useful down the road.

There is a similar relationship in continuous time relating the unit step function $u(t)$ as the integral of the Dirac $\delta(t)$.


How to perform model fitting for system identification

In any signal processing problem, there are usually two components: the signal model and the channel model.

The signal model is the mathematical description of how your ideal signal, call it $s(t)$ is generated.

The channel model is the mathematical description of how your channel corrupts or alters the signal.

One aim in telecommunications is to look at the received signal $r(t)$ --- $s(t)$ corrupted by the channel --- and try to undo the effects of the channel.

For example, suppose
$$
s(t) = A \sin(\omega_0 t + \phi)
$$
and suppose the channel is:
$$
r(t) = h(t)*s(t) + n(t)
$$
where $h(t)$ is the impulse response of the channel, and $n(t)$ is some additive noise.

In this context, your "AR"  filter is the $h(t)$.

If the context is in robotics or kinematics, then things get a little more complicated.

Did you open the "Algorithm" tab on the link in your comment? It's pretty straighforward:

  The polyfit MATLAB file forms the Vandermonde matrix, $\bf V$, whose elements are powers of $x$:
  $v_{i,j} = x_i^{n - j}$
  It then uses the backslash operator, $\backslash $, to solve the least squares problem ${\bf V} p \approx y$.


Does FFT on two samples from the same audio but with a very small sampling offset return the same results?

Fourier transforms possess a form of (in magnitude) shift-invariance, also called time-invariance, or shift property (see Properties of the Fourier Transform-shift properties): "a shift in time corresponds to a phase rotation in the frequency domain":
$$F\{x(t−t_0)\}=\exp(−j2πft_0)F\{x(t)\}\,.$$
This property of the continuous Fourier transform has analogues in its discretized versions. So many techniques to compensate for shifts or synchronize signals are based on this property. I'll name a few ones: cross-correlation, time-delay estimation, generalized cross-correlation (GCC), fractional delay filters. They can give you initial pointers.
In Matlab, functions like finddelay or alignsignals can be used.


How do I evaluate the performance of a segmentation algorithm?

Given that you are working only on Tanimoto Coefficient, i am trying to be more specific rather than giving generic answer with various different approaches.

The basic notation of Tanimoto Coefficient is as follows:

$$
T(A,B) = { N_{AB} \over N_A + N_B - N_{AB} }
$$

where $T$ is the desired result, over $A$ and $B$ images.

In this measure, we identify pixels either as belonging to a given segment i.e. it is segment pixel or it is background.
$N$ refers to number of pixels that are classified as segment pixel in the respective image. And $N_{AB}$ refers to the number of pixels which are classified as segment pixel in both images.

In this measure all pixels which qualifies as neither in A nor B is not calculated; only the pixels.

Also, both image must have same resolution and must have identical locations to for the segmented objects else even if the segmentation shape is right, the resultant overlap may not be right.

I am not getting into your MATLAB code, but here is the pseudo code that looks like.

Initialize N_A, N_B, N_AB;
for( all pixels @ x,y)
{
   if(image_A[x][y] == SEGMENT_CLASS_PIXEL)
    N_A += 1;

   if(image_B[x][y] == SEGMENT_CLASS_PIXEL)
    N_A += 1;

   if(image_A[x][y] == SEGMENT_CLASS_PIXEL
      && image_B[x][y] == SEGMENT_CLASS_PIXEL)
    N_AB += 1;

}

T = N_AB / (N_A + N_B - N_AB);


Implement the 3D Structure Tensor in C/C++?

There are many sides to this.

For typical kernel sizes, Fourier transform is not the most efficient method for convolution, if you are correctly utilizing the CPU for sliding window operations. That, of course, takes a good deal of software optimization. There are libraries to do that though (for example many of the deep learning libraries do 3D convolutions in CPU/GPU - I guess most of them do not use FFT). Given the specific symmetric nature of the Gaussian kernel, it is also very suitable for sliding operations. If the kernel size is very large, you could also try recursive convolutions.
The Fourier domain operations might be considered to be a better alternative if you have some luxury to pre-compute certain things. For example, if the image size is not changing and etc.
Probably the best answer would be to try and actually see it. I would prefer the separable convolutions. Vast amount of resources are also available for this.


Imposing sparse structure in Gabor space

In my experience, you do not necessarily need an orthogonal dictionary, see for example Compressed Sensing with Coherent and Redundant Dictionaries (DOI: 10.1016/j.acha.2010.10.002).

The paper Quantifying the performance of compressive
sensing on scalp EEG signals (DOI: 10.1109/ISABEL.2010.5702814) seems to be an example of the use of Gabor dictionaries (in EEG) in compressed sensing. I would just go ahead and try it.

If the chirp you are looking for can have a fairly fixed location in the Gabor domain, it might be possible to use compressive detection or classification instead of actual reconstruction. This is explained in Signal Processing With Compressive Measurements (10.1109/JSTSP.2009.2039178. You can see an example of how we have used it for symbol detection in IEEE 802.15.4 communication.


Consider a moving average over N samples- this is a simple FIR filter where each new output is the average of the past N samples. It is easy to see how high frequency noise can be filtered out (so is a low pass filter), and the longer time duration we include in the averaging window the lower will be the frequency cut off (just compare a stock market 30 day moving average to a 1 day moving average).

A moving average is a poor low pass filter, having a frequency response that approaches a Sinc function, which rolls off relatively slowly in frequency. By doing a weighted moving average where different samples are given different weights in the averaging process, we can significantly improve the frequency response - and coming up with the correct weights is the science of digital filter design.

IIR filters are similar except we are performing the average with previous outputs instead of past inputs.


3 dB cut-off frequency of first-order IIR high-pass filter

A discrete-time first-order high pass filter with unity gain at Nyquist and a zero at DC is described by the following difference equation:
$$y[n]=\frac{1+\alpha}{2}\big(x[n]-x[n-1]\big)+\alpha y[n-1],\qquad -1<\alpha<1\tag{1}$$
Its transfer function is given by
$$H(z)=\frac{1+\alpha}{2}\frac{1-z^{-1}}{1-\alpha z^{-1}}\tag{2}$$
Evaluating the squared magnitude of $(2)$ on the unit circle $z=e^{j\omega}$ and equating it to $\frac12$ ($-3$ dB) results in the following relation between $\alpha$ and the $3$ dB cut-off frequency $\omega_c$:
$$\begin{align}\big|H(e^{j\omega_c})\big|^2&=\frac{(1+\alpha)^2}{4}\frac{\left|1-e^{-j\omega_c}\right|^2}{\left|1-\alpha e^{-j\omega_c}\right|^2}\\&=\frac{(1+\alpha)^2}{4}\frac{2\big(1-\cos(\omega_c)\big)}{1-2\alpha\cos(\omega_c)+\alpha^2}\stackrel{!}{=}\frac12\tag{3}\end{align}$$
Eq. $(3)$ results in a quadratic equation for $\alpha$ with the solution
$$\alpha=\begin{cases}\displaystyle\frac{1-\sin(\omega_c)}{\cos(\omega_c)},&\omega_c\in(0,\pi)\setminus \frac{\pi}{2}\\0,&\omega_c=\frac{\pi}{2}\end{cases}\tag{4}$$
(where the requirement $|\alpha|<1$) has been taken into account).
For $\omega_c=\pi/2$ we obtain $\alpha=0$ and the corresponding filter is a simple $2$-tap FIR filter. All other cut-off frequencies $\omega_c\in(0,\pi)$ result in IIR filters.
The figure below shows the magnitude responses of $9$ high-pass filters with specified cut-off frequencies $0.1\pi,0.2\pi, \ldots,0.9\pi$. The corresponding values for $\alpha$ were computed according to Eq. $(4)$.


Output signal as convolution of impulse response and step signal

It is a normal IIR filter, and a very simple one at that.

Everything that can be described with a linear time invariant difference equation is an IIR filter and can be designed and analyzed using well established mathematical methods.

You calculate the output of the filter for ANY given input (impulse, step, the Dow Jones, Wagner's Nibelungen, etc.) with the equation
$$y[n] =\sum_{i=0}^{N}b_i\cdot x[n-i]-\sum_{i=1}^{N}a_i\cdot y[n-i]$$

In your case you have $b_0=0.2$ and $a_1 = -0.8$ with all the other a and b coefficients being zero, so the equation simplifies to
$$y[n] = 0.2 \cdot x[n] + 0.8 \cdot y[n-1] $$

That is indeed a simple first order low pass filter with a -3dB point at a normalized frequency of 0.0355 or thereabouts. It's not a very good low pass filter since it doesn't go to zero at the Nyquist frequency. A somewhat better choice would be

$$y[n] = 0.1 \cdot x[n] + 0.1 \cdot x[n-1] + 0.8 \cdot y[n-1] $$

There are plenty of smoothing algorithms that are not IIR filters but your specific example is just a simple first order low pass.

EDIT: based on the question in the comment
The equation above describes the operation for ALL input signals, not just an impulse.
Let's consider the input signal x = [1.0 1.1 1.05]; So we have x[0] = 1.0, x[1] = 1.1, x[2] = 1.05, and x[3...] = 0; You would calculate the output as follows

y[0] = 0.2*x[0] + 0.8*y[-1] = 0.2; // since y[-1] = 0;
y[1] = 0.2*x[1] + 0.8*y[0]  = 0.38;
y[2] = 0.2*x[2] + 0.8*y[1]  = 0.514;
y[3] =            0.8*y[2]  = 0.4112;     // since x[3] = 0;
y[4] =            0.8*y[3]  = 0.32896;
etc.


FFT peak estimation: Zero-padding vs signal repetitioon

If you repeat your data, you will have discontinuities where the sequence ends meet.  These will cause spurious peaks in your spectrum.  You could apply a window to the data to taper it at the edges, but that would be a fair amount of work, and you'd need to work out a window that doesn't color your spectrum.

I would just go with zero padding.  It will mess with the amplitude, but that should be fairly easy to deal with (scale by constant factor after FFT.)


Question related to Deconvolution and Fourier Transform

There could be 2 reasons why your 2 time series have different lengths :

1. They have a different sample rate

The situation

The sample rate of your signal is the frequency at which samples have been converted from the analog domain to the digital domain.
If 2 signals of the same duration have been sampled at different rates, the resulting time series will have different length( the longer being the one for the fastest sample rate).

The solution

The only solution to that problem is to resample one of your signals. For this, you will need an interpolation filter. In Matlab, the whole process is very simple. However, if you need to implement that in low-level code it's going to be more complex (lots of things to take in account). I won't go into technical details here, there's a lot you can find here and on the rest of the internet.

2. They have different durations

The situation

It's very much possible that your output signal has been recorded for a longer time than the input signal. In that case, you end up with 2 signals sampled at the same rate but with different lengths.

The solution

If that's the case, you're quite lucky as the solution is trivial. You just have to do something "zero-padding" which consists simply in adding as many zeroes at the end of the shortest signal as needed for it to have the same size at the other.

Additional info

There is however something that you have to consider with deconvolution : if your system has an infinite impulse response (IIR), your deconvolution will only give you an approached image of your transfer function. Think of it as for any function : if you only see part of the plot, the best you can do is to do educated guesses on the function's expression.

Moreover, depending on your sampling frequency and recording device, you might have a phenomenon of aliasing on your recorded response, and you will not see the frequency components that are above Nyquist's frequency.


I think you are seeking for an algorithm but this may help. If you don't know how $ s $ is distributed but you know its moments of first and second order, you can try with the linear estimator of the linear Minimum-mean Square Error (lmmse), with the following form,
$$
\hat{s}(y) = ay + b
$$
with
$$
a = \frac{ \mathrm{C}_{xy} }{ \sigma_{y}^{2} } \\
b = m_{s} - am_{y}
$$
where $ m $ denotes mean, $ \sigma $ variance and $ C $ covariance. It would help to know the moments of the noise, I mean, for example if is zero-mean and normalized to unity.


Why is the digital frequency highest at $\pi$ and what does it mean in terms of analog frequency?

Check this figure.

Consider the signal $x[n]=A$ which is constant but can be written as $x[n]=A\cos(0n)$, so its frequency is zero. There is no oscillation.

Let's increase the frequency. For the signal $x[n]=A\cos(\pi n/4)$, it oscillates at frequency $\omega_0=\pi/4$. The oscillation frequency increases.

For the frequency $\omega_0=\pi$, the signal is $x[n]=A\cos(\pi n)=A(-1)^n$, thus the signal oscillates between $A$ and $-A$ for each pair of consequtive samples. You can see that the signal cannot oscillate at a higher frequency, there is no faster oscillation than this.

For frequency $\omega_0=7\pi/4$, the signal can be written as $$x[n]=A\cos\Big(\frac{7\pi n}{4}\Big) = A\cos\Big(2\pi n-\frac{\pi n}{4}\Big)=A\cos\Big(\frac{\pi n}{4}\Big)$$ which oscillates at a lower frequency.

Finally, for $\omega_0=2\pi$, you have a constant signal again, no oscillation.

It is the nature of discrete time that introduces this effect.


Relation Between Compressed Sensing and Incoherence

The coherence can be used to bound the number of measurements $m$ required for successful reconstruction as seen for example in equation (6) in Candès & Wakin 2008:
$$m ≥ C \cdot \mu^2(\Phi, \Psi)\cdot S  \cdot \log n$$
where $C$ is some positive constant, $S$ is the sparsity, and $n$ is the size of the sparse signal.

One problem, however, is that this is not a very sharp bound. Donoho & Tanner 2010, for example, describe some brief history of this in Section X of the linked paper. I have not followed the development on this aspect of compressed sensing for a long time, so I do not know if there are newer or better explanations. So, in my opinion, the best you can use the coherence for in practice is as a rule of thumb: the less coherence, the better.


You're right, a causal and stable discrete-time system must have all its poles inside the unit circle. Since an allpass system with a zero $z_0$ inside the unit circle ($|z_0|<1$) must have a pole at $1/z_0^*$, i.e., outside the unit circle, it cannot be causal and stable.


Can somebody help me with some very basic questions about the fundamentals of information theory?

What is the nature of this relationship? There are many possibilities. Shannon used a logarithm because he wanted information to be additive, that is, $I(xy)=I(x)+I(y)$ for independent events $x$ and $y$.

You're correct that the arrival of a message at the receiver decreases its uncertainty. It may not be decreased to zero, because there may be an error in the communication. Also, a given channel can only support a certain information rate, which is its capacity.
It does not make sense to ask about messages of probability zero. The source alphabet is defined so that all symbols have probability greater than zero.


First of all, yes, what can be achieved by any machine learning tool is only as good as the dataset it was trained on.

The ML can't guess what properties of the high-dimensional data you want it to respect and which not.

So, you'll have to make sure all your data sources are either equal in their probability distributions when it comes to the feature you're interested in (language), essentially.

A relatively easy way of doing so is random sampling, where you categorize all data sources according to languages, and then for each language, pick a subset of recordings so that all sources are equally frequent.  You then train with different randomly picked subsets.

Really, this is Bayesian statistics basics; I don't think signal processing will help you here. You can sure still try to preprocess the audio data so that it only consists of relatively recording-method-invariate features (i.e. classical audio features, e.g. Cepstrum coefficients), but then you'd basically be doing what you'd instead wanted to use a neural net for: finding appropriate properties of the signal from which to tell the languages apart, and since the problem would persist, I'm pretty optimistic that the Resnet approach would still large correlate its classification to the data source, since some information would certainly persist through that feature-finding step, and your datasets still aren't even in terms of classification target pdf.


Cyclic prefix length in baseband OFDM with ZF equalizer

The channel which you have created is having 4-Taps and all taps are one after the other, meaning roughly there are 4 multi-paths and they are very close to each other. How close depends on what is the Sub-Carrier Spacing you would have assumed. Anyway, the point is only 3 samples of Cyclic Prefix would be enough and even with 0 or 2 samples of CP will not cause any noticeable degradation. The degradation is mostly due to AWGN which you are adding.

If you want to see the impact of 4-Tap channel, then I will suggest you to have larger delays multiplied to the magnitude of the four taps. One way to do it is as follows:

h = [.7*exp(1i*pi/4) zeros(1,50) .4*exp(1i*pi/3) zeros(1,75) .2*exp(1i*pi/6) zeros(1,100) .1*exp(1i*pi/5)]; % 4-taps channel

This will create four paths channel with 2nd, 3rd and 4rth path coming at 50, 126 and 227 samples respectively and hence the 2nd, 3rd and 4rth copy of one OFDM symbol will interfere with next OFDM symbol more.

With your current h assignment, all you will have is intra-symbol interference which can be delt with no CP since the channel is just 4-taps with no delays in between.

Also, consider removing AWGN and varying the amount of zeros I have added in between 4-Taps of h, to actually notice the impact of multipath inter-symbol interference and advantage of 512 length CP over 0/2 length CP.


How is the ends of a bandpass filtered signal reliable after applying the cosine taper?

A taper function is applied for spectral estimation if the artifacts from a non-tapered rectangular window provide less useful results that the artifacts due to the tapered window.

Not tapering a finite length sub-sample vector is equivalent to rectangular windowing, which convolves Sinc shaped artifacts onto the spectrum.


Not able to reach minimum phase using Hilbert transform

More taps. You don't have anywhere near enough taps for a filter that steep. Start large with 8192 or so cut to desired accuracy, if needed
Due to the low number of tabs you are seeing the effect of "circular" hilbert transform. See for example: http://andrewduncan.net/air/
How do you know you have zeros outside the unit circle? Calculating the roots of a high order polynomial is a numerically difficult problem and using a Matlab function like "roots()" will given you often inaccurate or plainly wrong results.


Audio Watermarking - What is typically embedded?

Research papers or textbooks won't delve into the specifics - they just assume you have a stream of binary data that needs to be reliably and undetectably embedded into the audio file.

What this bitstream is actually going to be is application dependent. It could be:

Digital rights management applications: A short serial number identifying the individual transaction in which the audio file has been sold. Given the small size (a hundred bits) of this information, it is typically repeated over and over the file to maximize the probability of recovering it.
Metadata management: Metadata such as artist name, song title, or even lyrics.
Proper steganography: Any secret data you want to share with the recipient of the audio file - without an eavesdroppers suspecting that you are exchanging anything other than music. This can be a text-file, an image - any binary data.
Informed source separation: Any information about the individual sources the audio file is made of, that would subsequently help in recovering these individual sources.
Second-screen interactive TV applications: These application allow a device such as a smartphone or a toy to display secondary information synchronized with a TV show. To allow the synchronization, the audio signal of the show contains a regular watermark containing a show identifier and time stamp.


matrix1 = rand(10, 100);
matrix2 = rand(10, 100);
for ind1=1:size(matrix1, 2)
    for ind2=1:size(matrix2, 2)
        [xc{ind1, ind2},lags{ind1, ind2}] = xcorr(matrix1(:,1),matrix2(:,2));
        % Here you can add any other operation you want
    end
end


As @jojek has already said in the comments, the formula reads
$$
P_\mathrm{x} = \dfrac{1}{T}\sum_{t=1}^{T}x^2(t)
$$
As $t$ usually denotes continous time I find this formulation a little bit odd. What your tutor probably actually meant is
$$
P_\mathrm{x} = \dfrac{1}{T}\int_{t=0}^{T}x^2(t) \mathrm{d}t
$$
It calculates the average power of the (continous-time) signal $x(t)$ of length $T$. In computer programs we're dealing with discrete-time (and discrete valued) signals, of course and the average power of the digital signal $x_n$ is given by
$$
\tilde P_\mathrm{x} = \dfrac{1}{N}\sum_{n=1}^{N}x_n^2,
$$
where $N$ is the length of $x_n$ in samples.


What's wrong with just doing sys = armax(data,[na,nb,nc,nk]); as per the Mathworks manual? sys is the variable. You could do sys[i] = ... to store several of them in an array (or perhaps sys{i}).

Then, to compare them, you could use getcov() and compare the covariances of each model.  These covariances tell you how "uncertain" the model is.

Another approach might be to look at the model residuals. If the residuals are white (uncorrelated), then the model is good.  The relatively energy / power of the residuals will be a measure of how accurate the model matches the data (the smaller the residual, the better the model).


It seems to be at a stadium which uses 12 cameras and Replay-Technologies special signal processing called FreeD: https://www.youtube.com/watch?v=N1kgt2VDjdM
The image you show seems to use 2-3 cameras, but the technology can work with multiple camera combinations from the 12 in 360'.
It's a particularly stunning real-time effect, and there are also very many 2D to 3D apps available for example in after effects and on google scholar.
Here is a list of replay-technologies videos:
https://www.youtube.com/results?search_query=replay+technologies+freed
https://www.tvtechnology.com/news/bringing-the-matrix-to-the-nfl


Ok so there are a couple ways to deal with this. As you pointed out there are things like the Costas Loop and PLL but I think band edge recovery is more intuitive. You can also continually apply STFT and track the frequency offset though you have the problem of time vs frequency localization.

I'm not an expert but I would probably square the signal and find the FFT peak, afterward I would do timing recovery and then apply decision directed feedback or a PLL after demodulation to track the remaining offset.

The band-edge filter timing recovery is pretty intuitive though fred harris derivation ignores some important things/obscures them for some reason. EX he tries to argue maximum-liklihood as a way to derive it but his derivation ignores issues like symbols not being time limited so it hides lots of stuff,


BPSK signal that has not been pulse shaped by any filter

To answer your question, NO!  A BPSK signal with a random data bit pattern that has not gone through any (additional) pulse shaping, meaning the pulse shape is rectangular in time equal to the bit duration, has a $Sinc^2$ power spectral density centered on the carrier frequency (as the Fourier Transform of a rectangular pulse is the Sinc function).

This signal has infinite bandwidth, but the power spectral density is dropping within an envelope of $1/f^2$ where f is the offset from the carrier (as expected from the $Sinc^2$ function). Any reduction of this bandwidth by filtering in the frequency domain, results in a finite transition time between symbols, so therefore induces a pulse shape in the time domain that is non-rectangular.

Another way to see this simply, a symbol that can transition from one state to the next in an infinitely short time (which is what the rectangular pulse shape is required to do) would require infinite bandwidth in order to transition that quickly. We slow down the transition rate on purpose in order to limit the bandwidth occupied since spectrum is such a valuable (and regulated) resource. We just do that very strategically to minimize inter-symbol interference, and optimize the SNR we can receive within the bandwidth used.


Filtering and Fourier Transforming, does the order matter?

There may be some slight differences due to the band-pass filtering of FFT results in the frequency domain being a circular convolution (with some wrap-around artifacts) rather than a pure linear convolution.

If you filter first, starting in time before your FFT window, any windowing artifacts from any out-of-band spectrum will be reduced before leaking into the FFT filter's pass-band.  Or alternatively, you could zero-pad the FFT window by the length of your filter's impulse response.


Why is the signum function $=2u(t) - 1$?

The signum function is defined by

$$\text{sgn}(t)=\begin{cases}-1,&t<0\\0,&t=0\\1,&t>0\end{cases}$$

Using the half-maximum convention, the unit step function is defined by

$$u(t)=\begin{cases}0,&t<0\\\frac12,&t=0\\1,&t>0\end{cases}$$

From these two definitions it should be obvious that

$$\text{sgn}(t)=2u(t)-1$$

must hold.